Study overview: This study was conducted at two location along the Po river in Northwestern Italy. The two sites were chosen to represent two distinct riparian condition. One site (Pian della Regina) was an alpine prairie with very sparse vegetation located above the treeline. The second site (Ostana) consisted of a dense riparian forest and much higher levels of shade and leaf inputs. Macroinvertebrate samples were hand collected and identified to species as well as classified to macroinvertebrate functional feeding groups.

After the surfaces of the macroinvertebrates were decontaminated, the internal bacterial communities were identified throught high throughtput sequencing of the 16S V4 region (515f, 806r) on an Illumina MiSeq platform. After filtering in QIIME2 (2018.11), bacterial reads were classified to taxanomic and functional orthologs (KEGG) using QIIME 2 and PiCRUSt 2.

A total of 26 invertebrates were sequenced (see metadata tab below for breakdown)

To see the corresponding code for each output click on the button labeled code above each box. All error bars are SEM and all multiple comparisons are adjusted with a FDR correction unless otherwise noted.

For code used to make figures see “ManuscriptFigures.r”. The .biom files and metadata are in the Github project repository.

Data overview

Sequencing

Based on the alpha rarefaction curves (Figure S1), samples were rarefied to 3,000 reads to limit the impact of differential library sizes on diversity metrics.

physeq
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 2420 taxa and 26 samples ]
## sample_data() Sample Data:       [ 26 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 2420 taxa by 14 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 2420 tips and 2400 internal nodes ]

Metadata

On average, predators, scrapers, and shredders were larger at the forested station (Ostana) than at the alpine prairie site.

Trtdata <- ddply(metadata, c("Sampling_station", "FFG","Taxon_name"), summarise,
                 N    = length(id),
                 meanWeight = mean(Mass),
                 sd   = sd(Mass),
                 se   = sd / sqrt(N)
)
Trtdata
##   Sampling_station      FFG           Taxon_name N meanWeight         sd
## 1           Ostana Filterer          Hydropsyche 4    41.4500  33.707615
## 2           Ostana Predator                Perla 2   190.5000 174.513954
## 3           Ostana  Scraper              Epeorus 5    22.7600   3.857849
## 4           Ostana Shredder Tipula_(Acutitipula) 3   534.9333 425.432889
## 5              PDR Predator          Dictyogenus 4    83.9500  16.270730
## 6              PDR  Scraper              Epeorus 4    11.6000   4.009988
## 7              PDR Shredder Tipula_(Acutitipula) 4    87.7000  42.817987
##           se
## 1  16.853808
## 2 123.400000
## 3   1.725283
## 4 245.623793
## 5   8.135365
## 6   2.004994
## 7  21.408993
ggplot(metadata,aes(x=Sampling_station,y=Mass,fill=Sampling_station))+geom_boxplot()+ #Write in labels from posthoc
  scale_fill_manual(values=cbPalette)+xlab("")+ylab("Mass (mg)")+labs(fill="FFG")+guides(fill=FALSE)+facet_wrap(~FFG,scales = "free_y")#+theme(axis.text.x = element_text(angle = 45, hjust = 1))+#+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters2), position=position_dodge(width=0.9), size=8,color="black")#+ 

head(metadata)
##     id FieldID           Taxon_name      FFG Sampling_station  Mass
## 1  II5    T1OS Tipula_(Acutitipula) Shredder           Ostana 989.1
## 2  II6    T2OS Tipula_(Acutitipula) Shredder           Ostana 470.0
## 3  II8    T4OS Tipula_(Acutitipula) Shredder           Ostana 145.7
## 4 II13    P1OS                Perla Predator           Ostana  67.1
## 5 II14    P2OS                Perla Predator           Ostana 313.9
## 6 II15    H2OS          Hydropsyche Filterer           Ostana   7.2
##    shannon faith_pd
## 1 7.086586 20.95169
## 2 7.430371 22.58281
## 3 7.383975 21.47860
## 4 6.952585 17.94331
## 5 6.200620 26.41718
## 6 6.153952 17.75811
metadataStatsSubset<- subset(metadata,FFG!="Predator")

metadataStatsSubset<- subset(metadataStatsSubset,FFG!="Filterer")
metadataStatsSubset
##      id FieldID           Taxon_name      FFG Sampling_station  Mass
## 1   II5    T1OS Tipula_(Acutitipula) Shredder           Ostana 989.1
## 2   II6    T2OS Tipula_(Acutitipula) Shredder           Ostana 470.0
## 3   II8    T4OS Tipula_(Acutitipula) Shredder           Ostana 145.7
## 10 E1OS    E1OS              Epeorus  Scraper           Ostana  20.2
## 11 E2OS    E2OS              Epeorus  Scraper           Ostana  27.5
## 12 E3OS    E3OS              Epeorus  Scraper           Ostana  19.0
## 13 E4OS    E4OS              Epeorus  Scraper           Ostana  20.8
## 14 E5OS    E5OS              Epeorus  Scraper           Ostana  26.3
## 15  II1     T1P Tipula_(Acutitipula) Shredder              PDR 131.8
## 16  II2     T2P Tipula_(Acutitipula) Shredder              PDR  53.9
## 17  II3     T3P Tipula_(Acutitipula) Shredder              PDR  48.2
## 18  II4     T4P Tipula_(Acutitipula) Shredder              PDR 116.9
## 23  E1P     E1P              Epeorus  Scraper              PDR  11.0
## 24  E3P     E3P              Epeorus  Scraper              PDR  17.4
## 25  E4P     E4P              Epeorus  Scraper              PDR   9.6
## 26  E5P     E5P              Epeorus  Scraper              PDR   8.4
##      shannon  faith_pd
## 1  7.0865864 20.951686
## 2  7.4303709 22.582812
## 3  7.3839746 21.478599
## 10 4.9093727 11.372657
## 11 1.9619653  9.905822
## 12 3.0828744  9.146140
## 13 4.8326128  7.641727
## 14 3.0611306  7.197870
## 15 7.2937707 21.728816
## 16 6.9992923 20.798101
## 17 6.0582633 11.170252
## 18 7.1744810 21.997409
## 23 1.6024907  4.266886
## 24 4.5283738  5.468116
## 25 0.7750815  4.605169
## 26 0.7526593  4.586619
print("Mass (mg)~ Sampling_station*FFG, method=ANOVA")
## [1] "Mass (mg)~ Sampling_station*FFG, method=ANOVA"
av=aov( Mass~ Sampling_station*FFG, data=metadataStatsSubset)
summary(av)
##                      Df Sum Sq Mean Sq F value  Pr(>F)   
## Sampling_station      1 109131  109131   3.563 0.08351 . 
## FFG                   1 319410  319410  10.427 0.00723 **
## Sampling_station:FFG  1 184026  184026   6.007 0.03054 * 
## Residuals            12 367594   30633                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print("(Since the interaction between sampling station and FFG was signinficant below are t-tests split out by FFG/n Mass (mg) ~Sampling-station (grouped by FFG, FDR correction)")
## [1] "(Since the interaction between sampling station and FFG was signinficant below are t-tests split out by FFG/n Mass (mg) ~Sampling-station (grouped by FFG, FDR correction)"
Scrapers<-subset(metadataStatsSubset, FFG=="Scraper")
t.test( Mass~ Sampling_station, data=Scrapers)
## 
##  Welch Two Sample t-test
## 
## data:  Mass by Sampling_station
## t = 4.2191, df = 6.4396, p-value = 0.004753
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##   4.79325 17.52675
## sample estimates:
## mean in group Ostana    mean in group PDR 
##                22.76                11.60
Shredders<-subset(metadataStatsSubset, FFG=="Shredder")
t.test( Mass~ Sampling_station, data=Shredders)
## 
##  Welch Two Sample t-test
## 
## data:  Mass by Sampling_station
## t = 1.8139, df = 2.0304, p-value = 0.2095
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -598.5184 1492.9851
## sample estimates:
## mean in group Ostana    mean in group PDR 
##             534.9333              87.7000
#compare_means(Mass_.mg. ~ Sampling_station,group.by = "FFG", data = metadata, p.adjust.method = "fdr",method="t.test")

Alpha Diversity

All three diversity metrics showed a similar pattern with scrapers having lower diversity than the other three groups.

Observed Species

plot_richness(physeq, x="FFG",color="FFG", measures=c("Observed"))+geom_boxplot(aes(x=FFG, y=value, color=FFG), alpha=0.05)+ylab("Observed Species")+geom_point(size=3)+guides(fill=FALSE)

Shannon Diversity

av=kruskal.test( shannon~ FFG, data=metadataStatsSubset) #fig.width = 8, fig.height = 8 for line above to change size
av2=kruskal.test( shannon~ Sampling_station, data=metadataStatsSubset) 
av
## 
##  Kruskal-Wallis rank sum test
## 
## data:  shannon by FFG
## Kruskal-Wallis chi-squared = 11.118, df = 1, p-value = 0.0008551
av2
## 
##  Kruskal-Wallis rank sum test
## 
## data:  shannon by Sampling_station
## Kruskal-Wallis chi-squared = 0.54044, df = 1, p-value = 0.4622
posthoc<-compare_means(shannon ~ FFG, data = metadataStatsSubset, p.adjust.method = "fdr")

Hyphenated<-as.character(paste0(posthoc$group1,"-",posthoc$group2))
difference<-posthoc$p.format
names(difference)<-Hyphenated

Letters<-multcompLetters(difference)

stats<-summarySE(metadataStatsSubset,measurevar="shannon",groupvars=c("FFG",na.rm=TRUE))
ggplot(stats,aes(x=FFG,y=shannon,fill=FFG))+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters), position=position_dodge(width=0.9), size=8,color="black")+ geom_bar(stat="identity")+  geom_errorbar(aes(ymin=shannon-se,ymax=shannon+se),color="black") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+ scale_fill_manual(values=cbPalette)+xlab("")+ylab("Shannon Diversity")+labs(fill="FFG")+guides(fill=FALSE)#+ theme(legend.justification=c(0.05,0.95), legend.position=c(0.05,0.95))

stats2<-summarySE(metadata,measurevar="shannon",groupvars=c("FFG","Sampling_station",na.rm=TRUE))
posthoc2<-compare_means(shannon ~ FFG, data = metadataStatsSubset, p.adjust.method = "fdr",group.by = "Sampling_station")


ggplot(metadata,aes(x=FFG,y=shannon,fill=FFG))+geom_boxplot() +#Write in labels from posthoc+geom_col()+geom_errorbar(aes(ymin=shannon-se,ymax=shannon+se),color="black")+
 scale_fill_manual(values=cbPalette)+xlab("")+ylab("Shannon Diversity")+labs(fill="FFG")+guides(fill=FALSE)+facet_wrap(~Sampling_station)#+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters2), position=position_dodge(width=0.9), size=8,color="black")#+ 

ggplot(metadata,aes(x=Sampling_station,y=shannon,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+xlab("")+ylab("Shannon Diversity")+labs(fill="FFG")+guides(fill=FALSE)+facet_wrap(~FFG)#+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters2), position=position_dodge(width=0.9), size=8,color="black")#+ 

av
## 
##  Kruskal-Wallis rank sum test
## 
## data:  shannon by FFG
## Kruskal-Wallis chi-squared = 11.118, df = 1, p-value = 0.0008551
posthoc
## # A tibble: 1 x 8
##   .y.     group1  group2         p p.adj p.format p.signif method  
##   <chr>   <chr>   <chr>      <dbl> <dbl> <chr>    <chr>    <chr>   
## 1 shannon Scraper Shredder 0.00103 0.001 0.001    **       Wilcoxon
Letters
##  Scraper Shredder 
##      "a"      "b"
stats
##        FFG na.rm N  shannon        sd        se        ci
## 1  Scraper    NA 9 2.834062 1.6633640 0.5544547 1.2785748
## 2 Shredder  TRUE 7 7.060963 0.4686418 0.1771300 0.4334214
posthoc2
## # A tibble: 2 x 9
##   Sampling_station .y.   group1 group2      p p.adj p.format p.signif
##   <fct>            <chr> <chr>  <chr>   <dbl> <dbl> <chr>    <chr>   
## 1 Ostana           shan~ Scrap~ Shred~ 0.0369 0.037 0.037    *       
## 2 PDR              shan~ Scrap~ Shred~ 0.0304 0.037 0.030    *       
## # ... with 1 more variable: method <chr>
stats2
##        FFG Sampling_station na.rm N  shannon        sd        se        ci
## 1 Filterer           Ostana    NA 4 6.502903 0.3062089 0.1531045 0.4872467
## 2 Predator           Ostana    NA 2 6.576603 0.5317195 0.3759825 4.7773101
## 3 Predator              PDR    NA 4 6.164048 1.3841309 0.6920654 2.2024611
## 4  Scraper           Ostana    NA 5 3.569591 1.2718220 0.5687761 1.5791756
## 5  Scraper              PDR    NA 4 1.914651 1.7867880 0.8933940 2.8431785
## 6 Shredder           Ostana  TRUE 3 7.300311 0.1865387 0.1076982 0.4633878
## 7 Shredder              PDR    NA 4 6.881452 0.5619605 0.2809802 0.8942045
av2
## 
##  Kruskal-Wallis rank sum test
## 
## data:  shannon by Sampling_station
## Kruskal-Wallis chi-squared = 0.54044, df = 1, p-value = 0.4622

AlphaDivSplitByFFG

print("Scrapers")
## [1] "Scrapers"
Scrapers<-subset(metadataStatsSubset, FFG=="Scraper")
shapiro.test(Scrapers$shannon)
## 
##  Shapiro-Wilk normality test
## 
## data:  Scrapers$shannon
## W = 0.8943, p-value = 0.2208
shapiro.test(Scrapers$faith_pd)
## 
##  Shapiro-Wilk normality test
## 
## data:  Scrapers$faith_pd
## W = 0.91387, p-value = 0.3439
t.test( shannon~ Sampling_station, data=Scrapers)
## 
##  Welch Two Sample t-test
## 
## data:  shannon by Sampling_station
## t = 1.5626, df = 5.2748, p-value = 0.1759
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -1.025406  4.335285
## sample estimates:
## mean in group Ostana    mean in group PDR 
##             3.569591             1.914651
t.test( faith_pd~ Sampling_station, data=Scrapers)
## 
##  Welch Two Sample t-test
## 
## data:  faith_pd by Sampling_station
## t = 5.3855, df = 4.8851, p-value = 0.003194
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  2.243935 6.398356
## sample estimates:
## mean in group Ostana    mean in group PDR 
##             9.052843             4.731698
kruskal.test( shannon~ Sampling_station, data=Scrapers)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  shannon by Sampling_station
## Kruskal-Wallis chi-squared = 2.94, df = 1, p-value = 0.08641
kruskal.test( faith_pd~ Sampling_station, data=Scrapers)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  faith_pd by Sampling_station
## Kruskal-Wallis chi-squared = 6, df = 1, p-value = 0.01431
print("Shredders")
## [1] "Shredders"
Shredders<-subset(metadataStatsSubset, FFG=="Shredder")
t.test(shannon~Sampling_station,data=Shredders)
## 
##  Welch Two Sample t-test
## 
## data:  shannon by Sampling_station
## t = 1.392, df = 3.8225, p-value = 0.2394
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.4321367  1.2698544
## sample estimates:
## mean in group Ostana    mean in group PDR 
##             7.300311             6.881452
t.test( faith_pd~ Sampling_station, data=Shredders)
## 
##  Welch Two Sample t-test
## 
## data:  faith_pd by Sampling_station
## t = 1.0402, df = 3.2033, p-value = 0.3703
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -5.364469 10.859245
## sample estimates:
## mean in group Ostana    mean in group PDR 
##             21.67103             18.92364
kruskal.test(shannon~Sampling_station,data=Shredders)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  shannon by Sampling_station
## Kruskal-Wallis chi-squared = 2, df = 1, p-value = 0.1573
kruskal.test( faith_pd~ Sampling_station, data=Shredders)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  faith_pd by Sampling_station
## Kruskal-Wallis chi-squared = 0.5, df = 1, p-value = 0.4795

AlphaDivSplitByStation

print("Ostana")
## [1] "Ostana"
Ostana<-subset(metadata, Sampling_station=="Ostana")
OstanaStatAlpha<-subset(Ostana,FFG!="Predator")
shapiro.test(OstanaStatAlpha$shannon)
## 
##  Shapiro-Wilk normality test
## 
## data:  OstanaStatAlpha$shannon
## W = 0.87409, p-value = 0.07366
shapiro.test(OstanaStatAlpha$faith_pd)
## 
##  Shapiro-Wilk normality test
## 
## data:  OstanaStatAlpha$faith_pd
## W = 0.8379, p-value = 0.02611
kruskal.test( shannon~ FFG, data=OstanaStatAlpha)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  shannon by FFG
## Kruskal-Wallis chi-squared = 9.6923, df = 2, p-value = 0.007859
compare_means(shannon ~ FFG, data = OstanaStatAlpha, p.adjust.method = "fdr")
## # A tibble: 3 x 8
##   .y.     group1   group2        p p.adj p.format p.signif method  
##   <chr>   <chr>    <chr>     <dbl> <dbl> <chr>    <chr>    <chr>   
## 1 shannon Filterer Scraper  0.0200 0.052 0.020    *        Wilcoxon
## 2 shannon Filterer Shredder 0.0518 0.052 0.052    ns       Wilcoxon
## 3 shannon Scraper  Shredder 0.0369 0.052 0.037    *        Wilcoxon
kruskal.test( faith_pd~ FFG, data=OstanaStatAlpha)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  faith_pd by FFG
## Kruskal-Wallis chi-squared = 8.2564, df = 2, p-value = 0.01611
compare_means(faith_pd ~ FFG, data = OstanaStatAlpha, p.adjust.method = "fdr")
## # A tibble: 3 x 8
##   .y.      group1   group2        p p.adj p.format p.signif method  
##   <chr>    <chr>    <chr>     <dbl> <dbl> <chr>    <chr>    <chr>   
## 1 faith_pd Filterer Scraper  0.0200 0.055 0.020    *        Wilcoxon
## 2 faith_pd Filterer Shredder 0.596  0.6   0.596    ns       Wilcoxon
## 3 faith_pd Scraper  Shredder 0.0369 0.055 0.037    *        Wilcoxon
print("PDR")
## [1] "PDR"
PDRStatAlpha<-subset(metadata, Sampling_station!="Ostana")
kruskal.test( shannon~ FFG, data=PDRStatAlpha)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  shannon by FFG
## Kruskal-Wallis chi-squared = 7.7308, df = 2, p-value = 0.02095
compare_means(shannon ~ FFG, data = PDRStatAlpha, p.adjust.method = "fdr")
## # A tibble: 3 x 8
##   .y.     group1   group2        p p.adj p.format p.signif method  
##   <chr>   <chr>    <chr>     <dbl> <dbl> <chr>    <chr>    <chr>   
## 1 shannon Predator Scraper  0.0304 0.046 0.03     *        Wilcoxon
## 2 shannon Predator Shredder 0.470  0.47  0.47     ns       Wilcoxon
## 3 shannon Scraper  Shredder 0.0304 0.046 0.03     *        Wilcoxon
kruskal.test( faith_pd~ FFG, data=PDRStatAlpha)
## 
##  Kruskal-Wallis rank sum test
## 
## data:  faith_pd by FFG
## Kruskal-Wallis chi-squared = 7.5385, df = 2, p-value = 0.02307
compare_means(faith_pd ~ FFG, data = PDRStatAlpha, p.adjust.method = "fdr")
## # A tibble: 3 x 8
##   .y.      group1   group2        p p.adj p.format p.signif method  
##   <chr>    <chr>    <chr>     <dbl> <dbl> <chr>    <chr>    <chr>   
## 1 faith_pd Predator Scraper  0.0304 0.046 0.03     *        Wilcoxon
## 2 faith_pd Predator Shredder 0.665  0.67  0.67     ns       Wilcoxon
## 3 faith_pd Scraper  Shredder 0.0304 0.046 0.03     *        Wilcoxon
posthoc<-compare_means(shannon ~ FFG, data = OstanaStatAlpha, p.adjust.method = "fdr")

Hyphenated<-as.character(paste0(posthoc$group1,"-",posthoc$group2))
difference<-posthoc$p.format
names(difference)<-Hyphenated

Letters<-multcompLetters(difference)

stats<-summarySE(OstanaStatAlpha,measurevar="shannon",groupvars=c("FFG",na.rm=TRUE))
ggplot(stats,aes(x=FFG,y=shannon,fill=FFG))+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters), position=position_dodge(width=0.9), size=8,color="black")+ geom_bar(stat="identity")+  geom_errorbar(aes(ymin=shannon-se,ymax=shannon+se),color="black") + scale_fill_manual(values=cbPalette)+xlab("")+ylab("Shannon Diversity (Ostana)")+labs(fill="FFG")+guides(fill=FALSE)#+ theme(legend.justification=c(0.05,0.95), legend.position=c(0.05,0.95))

posthoc<-compare_means(shannon ~ FFG, data = PDRStatAlpha, p.adjust.method = "fdr")

Hyphenated<-as.character(paste0(posthoc$group1,"-",posthoc$group2))
difference<-posthoc$p.format
names(difference)<-Hyphenated

Letters<-multcompLetters(difference)


stats<-summarySE(PDRStatAlpha,measurevar="shannon",groupvars=c("FFG",na.rm=TRUE))
ggplot(stats,aes(x=FFG,y=shannon,fill=FFG))+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters), position=position_dodge(width=0.9), size=8,color="black")+ geom_bar(stat="identity")+  geom_errorbar(aes(ymin=shannon-se,ymax=shannon+se),color="black") +  scale_fill_manual(values=cbPalette)+xlab("")+ylab("Shannon Diversity (PDR)")+labs(fill="FFG")+guides(fill=FALSE)#+ theme(legend.justification=c(0.05,0.95), legend.position=c(0.05,0.95))

posthoc<-compare_means(faith_pd ~ FFG, data = OstanaStatAlpha, p.adjust.method = "fdr")

Hyphenated<-as.character(paste0(posthoc$group1,"-",posthoc$group2))
difference<-posthoc$p.format
names(difference)<-Hyphenated

Letters<-multcompLetters(difference)

stats<-summarySE(OstanaStatAlpha,measurevar="faith_pd",groupvars=c("FFG",na.rm=TRUE))
ggplot(stats,aes(x=FFG,y=faith_pd,fill=FFG))+geom_text(aes(x=FFG, y=faith_pd+se+1,label=Letters$Letters), position=position_dodge(width=0.9), size=8,color="black")+ geom_bar(stat="identity")+  geom_errorbar(aes(ymin=faith_pd-se,ymax=faith_pd+se),color="black") +  scale_fill_manual(values=cbPalette)+xlab("")+ylab("Faith's PD (Ostana)")+labs(fill="FFG")+guides(fill=FALSE)#+ theme(legend.justification=c(0.05,0.95), legend.position=c(0.05,0.95))

posthoc<-compare_means(faith_pd ~ FFG, data = PDRStatAlpha, p.adjust.method = "fdr")

Hyphenated<-as.character(paste0(posthoc$group1,"-",posthoc$group2))
difference<-posthoc$p.format
names(difference)<-Hyphenated

Letters<-multcompLetters(difference)

stats<-summarySE(PDRStatAlpha,measurevar="faith_pd",groupvars=c("FFG",na.rm=TRUE))
ggplot(stats,aes(x=FFG,y=faith_pd,fill=FFG))+geom_text(aes(x=FFG, y=faith_pd+se+1,label=Letters$Letters), position=position_dodge(width=0.9), size=8,color="black")+ geom_bar(stat="identity")+  geom_errorbar(aes(ymin=faith_pd-se,ymax=faith_pd+se),color="black") +  scale_fill_manual(values=cbPalette)+xlab("")+ylab("Faith's PD (PDR)")+labs(fill="FFG")+guides(fill=FALSE)#+ theme(legend.justification=c(0.05,0.95), legend.position=c(0.05,0.95))

Faith’s PD

av=kruskal.test( faith_pd~ FFG, data=metadataStatsSubset) 
av2=kruskal.test( faith_pd~ Sampling_station, data=metadataStatsSubset) 
av
## 
##  Kruskal-Wallis rank sum test
## 
## data:  faith_pd by FFG
## Kruskal-Wallis chi-squared = 10.423, df = 1, p-value = 0.001245
av2
## 
##  Kruskal-Wallis rank sum test
## 
## data:  faith_pd by Sampling_station
## Kruskal-Wallis chi-squared = 0.89338, df = 1, p-value = 0.3446
posthoc<-compare_means(faith_pd ~ FFG, data = metadataStatsSubset, p.adjust.method = "fdr")


Hyphenated<-as.character(paste0(posthoc$group1,"-",posthoc$group2))
difference<-posthoc$p.format
names(difference)<-Hyphenated

Letters<-multcompLetters(difference)
stats<-summarySE(metadata,measurevar="faith_pd",groupvars=c("FFG",na.rm=TRUE))
ggplot(metadata,aes(x=FFG,y=faith_pd,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+xlab("")+ylab("Faith's PD")+labs(fill="FFG")#+ theme(legend.justification=c(0.05,0.95), legend.position=c(0.05,0.95))geom_text(aes(x=FFG, y=faith_pd+se+2,label=Letters$Letters), position=position_dodge(width=0.9), size=8,color="black")+ 

stats2<-summarySE(metadata,measurevar="faith_pd",groupvars=c("FFG","Sampling_station",na.rm=TRUE))
posthoc2<-compare_means(faith_pd ~ FFG, data = metadataStatsSubset, p.adjust.method = "fdr",group.by = "Sampling_station")


ggplot(metadata,aes(x=FFG,y=faith_pd,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+xlab("")+ylab("Faith PD")+labs(fill="FFG")+guides(fill=FALSE)+facet_wrap(~Sampling_station,scales = "free_x")+ theme(axis.text.x = element_text(angle = 45, hjust = 1))#+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters2), position=position_dodge(width=0.9), size=8,color="black")#+ 

ggplot(metadata,aes(x=Sampling_station,y=faith_pd,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+xlab("")+ylab("Faith's PD")+labs(fill="FFG")+guides(fill=FALSE)+facet_wrap(~FFG)#+geom_text(aes(x=FFG, y=shannon+se+1,label=Letters$Letters2), position=position_dodge(width=0.9), size=8,color="black")#+ 

av
## 
##  Kruskal-Wallis rank sum test
## 
## data:  faith_pd by FFG
## Kruskal-Wallis chi-squared = 10.423, df = 1, p-value = 0.001245
posthoc
## # A tibble: 1 x 8
##   .y.      group1  group2         p  p.adj p.format p.signif method  
##   <chr>    <chr>   <chr>      <dbl>  <dbl> <chr>    <chr>    <chr>   
## 1 faith_pd Scraper Shredder 0.00150 0.0015 0.0015   **       Wilcoxon
Letters
##  Scraper Shredder 
##      "a"      "b"
stats
##        FFG na.rm N  faith_pd       sd        se       ci
## 1 Filterer    NA 4 21.911376 2.894769 1.4473845 4.606223
## 2 Predator    NA 6 21.557631 5.684157 2.3205474 5.965157
## 3  Scraper    NA 9  7.132334 2.594210 0.8647367 1.994086
## 4 Shredder  TRUE 7 20.101096 3.984788 1.5061083 3.685314
posthoc2
## # A tibble: 2 x 9
##   Sampling_station .y.   group1 group2      p p.adj p.format p.signif
##   <fct>            <chr> <chr>  <chr>   <dbl> <dbl> <chr>    <chr>   
## 1 Ostana           fait~ Scrap~ Shred~ 0.0369 0.037 0.037    *       
## 2 PDR              fait~ Scrap~ Shred~ 0.0304 0.037 0.030    *       
## # ... with 1 more variable: method <chr>
av2
## 
##  Kruskal-Wallis rank sum test
## 
## data:  faith_pd by Sampling_station
## Kruskal-Wallis chi-squared = 0.89338, df = 1, p-value = 0.3446

Taxa Plots

GPr  = transform_sample_counts(physeq, function(x) x / sum(x) ) #transform samples based on relative abundance
GPr = filter_taxa(GPr, function(x) mean(x) > 1e-5, TRUE)
PhylumAll=tax_glom(GPr, "Phylum")

PhylumLevel = filter_taxa(PhylumAll, function(x) mean(x) > 1e-2, TRUE) #filter out any taxa lower tha 0.1%
FamilyAll=tax_glom(GPr,"Family")
FamilyLevel = filter_taxa(FamilyAll, function(x) mean(x) > 2e-2, TRUE) #filter out any taxa lower tha 1%
GenusAll=tax_glom(GPr,"Genus")
GenusLevel = filter_taxa(GenusAll, function(x) mean(x) > 2e-2, TRUE) #filter out any taxa lower tha 1%

Phylum Level

The top 7 phyla across all samples are given in the table below

df <- psmelt(PhylumLevel)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Phylum", "FFG"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
df2 <- psmelt(PhylumAll)
df2$Abundance<-df2$Abundance*100
PhylumAll<-ddply(df2, c("Phylum"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
PhylumSorted<-PhylumAll[order(-PhylumAll$mean),]
PhylumSorted[1:7,]
##             Phylum  N      mean        sd        se
## 20  Proteobacteria 26 51.896154 21.817321 4.2787287
## 5    Bacteroidetes 26 17.466667 14.180099 2.7809462
## 13      Firmicutes 26 13.275641 17.088823 3.3513939
## 23     Tenericutes 26 10.424359 19.936032 3.9097776
## 26 Verrucomicrobia 26  1.524359  1.846861 0.3621992
## 19  Planctomycetes 26  1.494872  1.825686 0.3580464
## 3   Actinobacteria 26  1.121795  1.100947 0.2159134
PhylumAll2<-ddply(df2, c("Phylum","FFG"), summarise, #For Rel abu tabs
                  N    = length(Abundance),
                  mean = mean(Abundance),
                  sd   = sd(Abundance),
                  se   = sd / sqrt(N)
)
#PhylumAll2 Run this to see all taxa, commented out since its a large list
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = Phylum),colour="black", stat="identity")+xlab("FFG")+ylab("Relative Abundance (> 1%)")+theme(axis.title.x=element_blank())+ theme(axis.text.x = element_text(angle = 45, hjust = 1))+ scale_fill_manual(values=cbPalette)
cdataplot#+  scale_x_discrete(labels=c("0 hrs", "24 hrs", "48 hrs","72 hrs"))

#position = "fill"#+geom_errorbar(aes(ymin=mean-se,ymax=mean+se),color="black",width=1)+geom_line(size=1.5,linetype="dashed")+geom_point(size=6)+ylab(
Trtdata2 <- ddply(df, c("Phylum", "FFG","Sampling_station"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
dfStats<-subset(df, FFG!="Predator")
dfStats<-subset(dfStats,FFG!="Filterer")
cdataplot2=ggplot(Trtdata2, aes(x=FFG,y=mean))+geom_bar(aes(fill = Phylum),colour="black", stat="identity")+xlab("FFG")+ylab("Relative Abundance (> 1%)")+theme(axis.title.x=element_blank())+ theme(axis.text.x = element_text(angle = 45, hjust = 1))+ scale_fill_manual(values=cbPalette)+facet_wrap(~Sampling_station,scales = "free_x")
cdataplot2

########
compare_means(Abundance ~ FFG, data = df, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test")
## # A tibble: 7 x 7
##   Phylum          .y.              p  p.adj p.format p.signif method       
##   <fct>           <chr>        <dbl>  <dbl> <chr>    <chr>    <chr>        
## 1 Proteobacteria  Abundance 0.0147   0.021  0.01466  *        Kruskal-Wall~
## 2 Tenericutes     Abundance 0.0643   0.075  0.06432  ns       Kruskal-Wall~
## 3 Bacteroidetes   Abundance 0.000400 0.0014 0.00040  ***      Kruskal-Wall~
## 4 Firmicutes      Abundance 0.000209 0.0014 0.00021  ***      Kruskal-Wall~
## 5 Planctomycetes  Abundance 0.0127   0.021  0.01274  *        Kruskal-Wall~
## 6 Verrucomicrobia Abundance 0.00184  0.0043 0.00184  **       Kruskal-Wall~
## 7 Actinobacteria  Abundance 0.118    0.12   0.11762  ns       Kruskal-Wall~
Means<-compare_means(Abundance ~ FFG, data = df, group.by = "Phylum", p.adjust.method = "fdr")


Trtdata <- ddply(df, c("Phylum", "FFG"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
SigList<-length(unique(Trtdata$Phylum))
for (i in levels(Means$Phylum)){
  Tax<-i
  TaxAbundance<-subset(Means,Phylum==i )
  Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
  difference<-TaxAbundance$p
  names(difference)<-Hyphenated
  Letters<-multcompLetters(difference)
  #print(Letters)
  SigList[i]<-Letters
  
}
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
vec<-unlist(SigList)
vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Family)+xlab("FFG")+ylab("Relative Abundance (> 1%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Phylum)+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+geom_text(aes(x=FFG, y=mean+se+10,label=vec))  +scale_fill_manual(values=cbPalette)
cdataplot

head(Trtdata)
##           Phylum      FFG N       mean         sd        se
## 1 Actinobacteria Filterer 4  0.3833333  0.4811252 0.2405626
## 2 Actinobacteria Predator 6  1.3111111  0.3862450 0.1576839
## 3 Actinobacteria  Scraper 9  1.2259259  1.2891403 0.4297134
## 4 Actinobacteria Shredder 7  1.2476190  1.4698018 0.5555329
## 5  Bacteroidetes Filterer 4 22.6166667  3.7669616 1.8834808
## 6  Bacteroidetes Predator 6 32.7000000 10.7146006 4.3742174
dfStatsPlot<-subset(Trtdata,Phylum!="Actinobacteria"&Phylum!= "Verrucomicrobia")
#vec<-c("a","b","a","b","a","b","a","b","a","b")
dfStatsPlot # Only show significant taxa 
##            Phylum      FFG N       mean         sd         se
## 5   Bacteroidetes Filterer 4 22.6166667  3.7669616 1.88348082
## 6   Bacteroidetes Predator 6 32.7000000 10.7146006 4.37421739
## 7   Bacteroidetes  Scraper 9  0.8444444  0.6502136 0.21673788
## 8   Bacteroidetes Shredder 7 22.8380952  6.0532033 2.28789578
## 9      Firmicutes Filterer 4 20.0250000 13.0079285 6.50396425
## 10     Firmicutes Predator 6  0.2555556  0.4385160 0.17902341
## 11     Firmicutes  Scraper 9  0.7888889  1.2328828 0.41096093
## 12     Firmicutes Shredder 7 36.6333333  7.7272486 2.92062543
## 13 Planctomycetes Filterer 4  1.3083333  0.4085884 0.20429418
## 14 Planctomycetes Predator 6  1.6000000  0.7636171 0.31174539
## 15 Planctomycetes  Scraper 9  2.4444444  2.7655118 0.92183727
## 16 Planctomycetes Shredder 7  0.2904762  0.2052228 0.07756693
## 17 Proteobacteria Filterer 4 45.1250000 15.0889554 7.54447768
## 18 Proteobacteria Predator 6 51.1833333  9.3604072 3.82137022
## 19 Proteobacteria  Scraper 9 68.0037037 27.1006890 9.03356301
## 20 Proteobacteria Shredder 7 35.6666667  9.4683488 3.57869948
## 21    Tenericutes Filterer 4  3.0750000  3.1347958 1.56739788
## 22    Tenericutes Predator 6  5.0666667 10.4961580 4.28503857
## 23    Tenericutes  Scraper 9 24.9888889 28.2112665 9.40375549
## 24    Tenericutes Shredder 7  0.4904762  0.3207135 0.12121831
cdataplot=ggplot(dfStatsPlot, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Family)+xlab("FFG")+ylab("Relative Abundance (> 1%)") + theme(axis.text.x = element_text(angle = 0, hjust = 0.5))+theme(axis.title.x=element_blank())+facet_wrap(~Phylum)+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+#geom_text(aes(x=FFG, y=mean+se+10,label=vec))+
  scale_fill_manual(values=cbPalette)+theme(legend.justification=c(0.05,0.95), legend.position=c(0.7,0.4))
cdataplot

#NComparisons<-length(unique(metadata$FFG))*length(unique(Trtdata$Phylum))
#SigList<-length(unique(Trtdata$Phylum))
#SigLetters2<-vector(length=NComparisons)
#vec<-unlist(lst)
Means=compare_means(Abundance ~ FFG, data = dfStats, group.by = "Phylum", p.adjust.method = "fdr")
MeansStation=compare_means(Abundance ~ Sampling_station, data = dfStats, group.by = "Phylum", p.adjust.method = "fdr")

Means
## # A tibble: 7 x 9
##   Phylum     .y.     group1  group2       p  p.adj p.format p.signif method
##   <fct>      <chr>   <chr>   <chr>    <dbl>  <dbl> <chr>    <chr>    <chr> 
## 1 Proteobac~ Abunda~ Scraper Shred~ 0.0149  0.021  0.015    *        Wilco~
## 2 Tenericut~ Abunda~ Scraper Shred~ 0.0148  0.021  0.015    *        Wilco~
## 3 Firmicutes Abunda~ Scraper Shred~ 0.00103 0.0036 0.001    **       Wilco~
## 4 Bacteroid~ Abunda~ Scraper Shred~ 0.00103 0.0036 0.001    **       Wilco~
## 5 Planctomy~ Abunda~ Scraper Shred~ 0.0148  0.021  0.015    *        Wilco~
## 6 Actinobac~ Abunda~ Scraper Shred~ 0.751   0.75   0.751    ns       Wilco~
## 7 Verrucomi~ Abunda~ Scraper Shred~ 0.289   0.34   0.289    ns       Wilco~
MeansStation
## # A tibble: 7 x 9
##   Phylum       .y.     group1 group2      p p.adj p.format p.signif method 
##   <fct>        <chr>   <chr>  <chr>   <dbl> <dbl> <chr>    <chr>    <chr>  
## 1 Proteobacte~ Abunda~ Ostana PDR    0.958   1    0.958    ns       Wilcox~
## 2 Tenericutes  Abunda~ Ostana PDR    0.528   0.89 0.528    ns       Wilcox~
## 3 Firmicutes   Abunda~ Ostana PDR    1       1    1.000    ns       Wilcox~
## 4 Bacteroidet~ Abunda~ Ostana PDR    0.637   0.89 0.637    ns       Wilcox~
## 5 Planctomyce~ Abunda~ Ostana PDR    0.0180  0.13 0.018    *        Wilcox~
## 6 Actinobacte~ Abunda~ Ostana PDR    0.495   0.89 0.495    ns       Wilcox~
## 7 Verrucomicr~ Abunda~ Ostana PDR    0.0738  0.26 0.074    ns       Wilcox~
Means<-compare_means(Abundance ~ FFG, data = dfStats, group.by = "Phylum", p.adjust.method = "fdr")


Trtdata <- ddply(dfStats, c("Phylum", "FFG"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
SigList<-length(unique(Trtdata$Phylum))
for (i in levels(Means$Phylum)){
  Tax<-i
  TaxAbundance<-subset(Means,Phylum==i )
  Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
  difference<-TaxAbundance$p
  names(difference)<-Hyphenated
  Letters<-multcompLetters(difference)
  #print(Letters)
  SigList[i]<-Letters
  
}
vec<-unlist(SigList)
vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Family)+xlab("FFG")+ylab("Relative Abundance (> 1%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Phylum)+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+geom_text(aes(x=FFG, y=mean+se+5,label=vec))+scale_fill_manual(values=cbPalette)
cdataplot

#######

Split by Station

Since there were not enough samples collect to compare between all groups at both stations, only the groups present with greater than 3 samples were compared at each station.

Ostana (Forest)

df <- psmelt(PhylumLevel)
df$Abundance=df$Abundance*100
dfOstana<-subset(df,Sampling_station =="Ostana")
dfStatsPDR<-subset(df,Sampling_station!="Ostana")


dfStatsOstana<-subset(dfOstana,FFG!="Predator")

compare_means(Abundance ~ FFG, data = dfStatsOstana, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test")
## # A tibble: 7 x 7
##   Phylum          .y.             p p.adj p.format p.signif method        
##   <fct>           <chr>       <dbl> <dbl> <chr>    <chr>    <chr>         
## 1 Proteobacteria  Abundance 0.291   0.290 0.2906   ns       Kruskal-Wallis
## 2 Tenericutes     Abundance 0.0549  0.077 0.0549   ns       Kruskal-Wallis
## 3 Firmicutes      Abundance 0.0121  0.042 0.0121   *        Kruskal-Wallis
## 4 Bacteroidetes   Abundance 0.00990 0.042 0.0099   **       Kruskal-Wallis
## 5 Planctomycetes  Abundance 0.0430  0.077 0.0430   *        Kruskal-Wallis
## 6 Actinobacteria  Abundance 0.109   0.13  0.1089   ns       Kruskal-Wallis
## 7 Verrucomicrobia Abundance 0.0484  0.077 0.0484   *        Kruskal-Wallis
Means<-compare_means(Abundance ~ FFG, data = dfStatsOstana, group.by = "Phylum", p.adjust.method = "fdr")

keeps <- c("Phylum","group1","group2","p.adj","method","p.signif")
keeps=Means[keeps]

test3 <- list('Phylum'= keeps$Phylum,'group1'=keeps$group1,'group2'= keeps$group2,'p.adj'=keeps$p.adj)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p.adj>0.1),]            
FilteredResults
## [1] Phylum group1 group2 p.adj 
## <0 rows> (or 0-length row.names)

Pian della Regina (Alpine Prairie)

print("Pian della Regina")
## [1] "Pian della Regina"
Means<-compare_means(Abundance ~ FFG, data = dfStatsPDR, group.by = "Phylum", p.adjust.method = "fdr")

keeps <- c("Phylum","group1","group2","p.adj","method","p.signif")
keeps=Means[keeps]

test3 <- list('Phylum'= keeps$Phylum,'group1'=keeps$group1,'group2'= keeps$group2,'p.adj'=keeps$p.adj)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p.adj>0.1),]            
FilteredResults
##             Phylum   group1   group2 p.adj
## 2   Proteobacteria Predator Shredder 0.071
## 3   Proteobacteria  Scraper Shredder 0.071
## 7    Bacteroidetes Predator  Scraper 0.071
## 9    Bacteroidetes  Scraper Shredder 0.071
## 11      Firmicutes Predator Shredder 0.071
## 12      Firmicutes  Scraper Shredder 0.071
## 13 Verrucomicrobia Predator  Scraper 0.071
## 14 Verrucomicrobia Predator Shredder 0.071
## 20  Planctomycetes Predator Shredder 0.071
#########

############
ggplot(dfOstana, aes(x=FFG,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Phylum)+ylab("Relative Abundance Ostana(>1%)")

ggplot(dfStatsPDR, aes(x=FFG,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Phylum)+ylab("Relative Abundance PDR(>1%)")

print("Family level")
## [1] "Family level"
df <- psmelt(FamilyLevel)
df$Abundance=df$Abundance*100
dfOstana<-subset(df,Sampling_station =="Ostana")
dfStatsPDR<-subset(df,Sampling_station!="Ostana")

dfStatsOstana<-subset(dfOstana,FFG!="Predator")

print("Ostana")
## [1] "Ostana"
compare_means(Abundance ~ FFG, data = dfStatsOstana, group.by = "Family", p.adjust.method = "fdr",method="kruskal.test")
## # A tibble: 10 x 7
##    Family            .y.             p p.adj p.format p.signif method      
##    <fct>             <chr>       <dbl> <dbl> <chr>    <chr>    <chr>       
##  1 Aeromonadaceae    Abundance 0.0817  0.091 0.0817   ns       Kruskal-Wal~
##  2 Comamonadaceae    Abundance 0.00786 0.021 0.0079   **       Kruskal-Wal~
##  3 Rhodobacteraceae  Abundance 0.886   0.89  0.8865   ns       Kruskal-Wal~
##  4 Ruminococcaceae   Abundance 0.00712 0.021 0.0071   **       Kruskal-Wal~
##  5 Desulfovibrionac~ Abundance 0.0105  0.021 0.0105   *        Kruskal-Wal~
##  6 Lachnospiraceae   Abundance 0.0118  0.021 0.0118   *        Kruskal-Wal~
##  7 Rikenellaceae     Abundance 0.0330  0.041 0.0330   *        Kruskal-Wal~
##  8 Methylophilaceae  Abundance 0.0139  0.021 0.0139   *        Kruskal-Wal~
##  9 Enterobacteriace~ Abundance 0.0144  0.021 0.0144   *        Kruskal-Wal~
## 10 Cytophagaceae     Abundance 0.0139  0.021 0.0139   *        Kruskal-Wal~
Means<-compare_means(Abundance ~ FFG, data = dfStatsOstana, group.by = "Family", p.adjust.method = "fdr")

keeps <- c("Family","group1","group2","p.adj","method","p.signif")
keeps=Means[keeps]

test3 <- list('Family'= keeps$Family,'group1'=keeps$group1,'group2'= keeps$group2,'p.adj'=keeps$p.adj)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p.adj>0.1),]            
FilteredResults
##                 Family   group1   group2 p.adj
## 3       Comamonadaceae Filterer  Scraper 0.058
## 4       Comamonadaceae Filterer Shredder 0.088
## 5       Comamonadaceae  Scraper Shredder 0.076
## 9      Ruminococcaceae Filterer  Scraper 0.054
## 11     Ruminococcaceae  Scraper Shredder 0.054
## 12 Desulfovibrionaceae Filterer  Scraper 0.054
## 14 Desulfovibrionaceae  Scraper Shredder 0.054
## 15     Lachnospiraceae Filterer  Scraper 0.054
## 17     Lachnospiraceae  Scraper Shredder 0.054
## 18       Rikenellaceae Filterer  Scraper 0.079
## 20       Rikenellaceae  Scraper Shredder 0.054
## 21    Methylophilaceae Filterer  Scraper 0.054
## 23    Methylophilaceae  Scraper Shredder 0.063
## 25  Enterobacteriaceae Filterer Shredder 0.079
## 26  Enterobacteriaceae  Scraper Shredder 0.076
## 27       Cytophagaceae Filterer  Scraper 0.054
## 29       Cytophagaceae  Scraper Shredder 0.063
print("Pian della Regina")
## [1] "Pian della Regina"
compare_means(Abundance ~ FFG, data = dfStatsPDR, group.by = "Family", p.adjust.method = "fdr",method="kruskal.test")
## # A tibble: 10 x 7
##    Family            .y.             p p.adj p.format p.signif method      
##    <fct>             <chr>       <dbl> <dbl> <chr>    <chr>    <chr>       
##  1 Enterobacteriace~ Abundance 0.0366  0.041 0.0366   *        Kruskal-Wal~
##  2 Cytophagaceae     Abundance 0.00609 0.013 0.0061   **       Kruskal-Wal~
##  3 Desulfovibrionac~ Abundance 0.00537 0.013 0.0054   **       Kruskal-Wal~
##  4 Lachnospiraceae   Abundance 0.00921 0.013 0.0092   **       Kruskal-Wal~
##  5 Rikenellaceae     Abundance 0.00537 0.013 0.0054   **       Kruskal-Wal~
##  6 Rhodobacteraceae  Abundance 0.0249  0.031 0.0249   *        Kruskal-Wal~
##  7 Comamonadaceae    Abundance 0.00715 0.013 0.0072   **       Kruskal-Wal~
##  8 Ruminococcaceae   Abundance 0.00921 0.013 0.0092   **       Kruskal-Wal~
##  9 Aeromonadaceae    Abundance 0.368   0.37  0.3679   ns       Kruskal-Wal~
## 10 Methylophilaceae  Abundance 0.00609 0.013 0.0061   **       Kruskal-Wal~
Means<-compare_means(Abundance ~ FFG, data = dfStatsPDR, group.by = "Family", p.adjust.method = "fdr")

keeps <- c("Family","group1","group2","p.adj","method","p.signif")
keeps=Means[keeps]

test3 <- list('Family'= keeps$Family,'group1'=keeps$group1,'group2'= keeps$group2,'p.adj'=keeps$p.adj)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p.adj>0.1),]            
FilteredResults
##                 Family   group1   group2 p.adj
## 3   Enterobacteriaceae  Scraper Shredder 0.041
## 4        Cytophagaceae Predator  Scraper 0.041
## 5        Cytophagaceae Predator Shredder 0.041
## 6        Cytophagaceae  Scraper Shredder 0.041
## 7  Desulfovibrionaceae Predator Shredder 0.041
## 8  Desulfovibrionaceae  Scraper Shredder 0.041
## 10     Lachnospiraceae Predator Shredder 0.041
## 11     Lachnospiraceae  Scraper Shredder 0.041
## 12       Rikenellaceae Predator Shredder 0.041
## 13       Rikenellaceae  Scraper Shredder 0.041
## 14    Rhodobacteraceae Predator  Scraper 0.041
## 15    Rhodobacteraceae Predator Shredder 0.041
## 17      Comamonadaceae Predator  Scraper 0.041
## 18      Comamonadaceae Predator Shredder 0.041
## 19      Comamonadaceae  Scraper Shredder 0.041
## 21     Ruminococcaceae Predator Shredder 0.041
## 22     Ruminococcaceae  Scraper Shredder 0.041
## 25    Methylophilaceae Predator  Scraper 0.041
## 26    Methylophilaceae Predator Shredder 0.041
## 27    Methylophilaceae  Scraper Shredder 0.041
ggplot(dfOstana, aes(x=FFG,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Family)+ylab("Relative Abundance Ostana(>2%)")

ggplot(dfStatsPDR, aes(x=FFG,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Family)+ylab("Relative Abundance PDR(>2%)")

print("Genus Level")
## [1] "Genus Level"
df <- psmelt(GenusLevel)
df$Abundance=df$Abundance*100
dfOstana<-subset(df,Sampling_station =="Ostana")
dfStatsPDR<-subset(df,Sampling_station!="Ostana")


dfStatsOstana<-subset(dfOstana,FFG!="Predator")

print("Ostana")
## [1] "Ostana"
compare_means(Abundance ~ FFG, data = dfStatsOstana, group.by = "Genus", p.adjust.method = "fdr",method="kruskal.test")
## # A tibble: 3 x 7
##   Genus         .y.            p p.adj p.format p.signif method        
##   <fct>         <chr>      <dbl> <dbl> <chr>    <chr>    <chr>         
## 1 Rhodobacter   Abundance 0.886  0.89  0.886    ns       Kruskal-Wallis
## 2 Clostridium   Abundance 0.0118 0.021 0.012    *        Kruskal-Wallis
## 3 Methylotenera Abundance 0.0139 0.021 0.014    *        Kruskal-Wallis
compare_means(Abundance ~ FFG, data = dfStatsOstana, group.by = "Genus", p.adjust.method = "fdr")
## # A tibble: 9 x 9
##   Genus      .y.     group1  group2       p p.adj p.format p.signif method 
##   <fct>      <chr>   <chr>   <chr>    <dbl> <dbl> <chr>    <chr>    <chr>  
## 1 Rhodobact~ Abunda~ Filter~ Scraper 0.903  0.9   0.903    ns       Wilcox~
## 2 Rhodobact~ Abunda~ Filter~ Shredd~ 0.860  0.9   0.860    ns       Wilcox~
## 3 Rhodobact~ Abunda~ Scraper Shredd~ 0.766  0.9   0.766    ns       Wilcox~
## 4 Clostridi~ Abunda~ Filter~ Scraper 0.0108 0.05  0.011    *        Wilcox~
## 5 Clostridi~ Abunda~ Filter~ Shredd~ 0.596  0.89  0.596    ns       Wilcox~
## 6 Clostridi~ Abunda~ Scraper Shredd~ 0.0168 0.05  0.017    *        Wilcox~
## 7 Methylote~ Abunda~ Filter~ Scraper 0.0151 0.05  0.015    *        Wilcox~
## 8 Methylote~ Abunda~ Filter~ Shredd~ 0.596  0.89  0.596    ns       Wilcox~
## 9 Methylote~ Abunda~ Scraper Shredd~ 0.0262 0.059 0.026    *        Wilcox~
print("Pian della Regina")
## [1] "Pian della Regina"
compare_means(Abundance ~ FFG, data = dfStatsPDR, group.by = "Genus", p.adjust.method = "fdr",method="kruskal.test")
## # A tibble: 3 x 7
##   Genus         .y.             p  p.adj p.format p.signif method        
##   <fct>         <chr>       <dbl>  <dbl> <chr>    <chr>    <chr>         
## 1 Clostridium   Abundance 0.00537 0.0091 0.0054   **       Kruskal-Wallis
## 2 Rhodobacter   Abundance 0.0246  0.025  0.0246   *        Kruskal-Wallis
## 3 Methylotenera Abundance 0.00609 0.0091 0.0061   **       Kruskal-Wallis
compare_means(Abundance ~ FFG, data = dfStatsPDR, group.by = "Genus", p.adjust.method = "fdr")
## # A tibble: 8 x 9
##   Genus      .y.     group1  group2       p p.adj p.format p.signif method 
##   <fct>      <chr>   <chr>   <chr>    <dbl> <dbl> <chr>    <chr>    <chr>  
## 1 Clostridi~ Abunda~ Predat~ Shredd~ 0.0211 0.035 0.021    *        Wilcox~
## 2 Clostridi~ Abunda~ Scraper Shredd~ 0.0211 0.035 0.021    *        Wilcox~
## 3 Rhodobact~ Abunda~ Predat~ Scraper 0.0304 0.035 0.030    *        Wilcox~
## 4 Rhodobact~ Abunda~ Predat~ Shredd~ 0.0294 0.035 0.029    *        Wilcox~
## 5 Rhodobact~ Abunda~ Scraper Shredd~ 1      1     1.000    ns       Wilcox~
## 6 Methylote~ Abunda~ Predat~ Scraper 0.0211 0.035 0.021    *        Wilcox~
## 7 Methylote~ Abunda~ Predat~ Shredd~ 0.0304 0.035 0.030    *        Wilcox~
## 8 Methylote~ Abunda~ Scraper Shredd~ 0.0211 0.035 0.021    *        Wilcox~
ggplot(dfOstana, aes(x=FFG,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Genus)+ylab("Relative Abundance Ostana(>2%)")

ggplot(dfStatsPDR, aes(x=FFG,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Genus)+ylab("Relative Abundance PDR(>2%)")

Split by FFG

This set of tests is comparing bacterial between sites for a particular group (e.g. Shredders at Ostana va PDR)

df <- psmelt(PhylumLevel)
df$Abundance=df$Abundance*100
dfShredders<-subset(df,FFG =="Shredder")
dfScrapers<-subset(df,FFG == "Scraper")

print("Shredders")
## [1] "Shredders"
compare_means(Abundance ~ Sampling_station, data = dfShredders, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test")
## # A tibble: 7 x 7
##   Phylum          .y.            p p.adj p.format p.signif method        
##   <fct>           <chr>      <dbl> <dbl> <chr>    <chr>    <chr>         
## 1 Proteobacteria  Abundance 0.480   0.56 0.480    ns       Kruskal-Wallis
## 2 Firmicutes      Abundance 0.480   0.56 0.480    ns       Kruskal-Wallis
## 3 Bacteroidetes   Abundance 0.0339  0.18 0.034    *        Kruskal-Wallis
## 4 Actinobacteria  Abundance 0.0771  0.18 0.077    ns       Kruskal-Wallis
## 5 Tenericutes     Abundance 0.372   0.56 0.372    ns       Kruskal-Wallis
## 6 Planctomycetes  Abundance 0.0771  0.18 0.077    ns       Kruskal-Wallis
## 7 Verrucomicrobia Abundance 0.589   0.59 0.589    ns       Kruskal-Wallis
Means<-compare_means(Abundance ~ Sampling_station, data = dfShredders, group.by = "Phylum", p.adjust.method = "fdr")

keeps <- c("Phylum","group1","group2","p","method","p.signif")
keeps=Means[keeps]

test3 <- list('Phylum'= keeps$Phylum,'group1'=keeps$group1,'group2'= keeps$group2,'p'=keeps$p)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p>0.1),]            
FilteredResults
##          Phylum group1 group2          p
## 3 Bacteroidetes Ostana    PDR 0.05182993
print("Scrapers")
## [1] "Scrapers"
compare_means(Abundance ~ Sampling_station, data = dfScrapers, group.by = "Phylum", p.adjust.method = "fdr",method="kruskal.test")
## # A tibble: 7 x 7
##   Phylum          .y.            p p.adj p.format p.signif method        
##   <fct>           <chr>      <dbl> <dbl> <chr>    <chr>    <chr>         
## 1 Proteobacteria  Abundance 0.327  0.570 0.33     ns       Kruskal-Wallis
## 2 Tenericutes     Abundance 0.462  0.62  0.46     ns       Kruskal-Wallis
## 3 Planctomycetes  Abundance 0.0500 0.35  0.05     ns       Kruskal-Wallis
## 4 Actinobacteria  Abundance 0.623  0.62  0.62     ns       Kruskal-Wallis
## 5 Firmicutes      Abundance 0.623  0.62  0.62     ns       Kruskal-Wallis
## 6 Bacteroidetes   Abundance 0.221  0.51  0.22     ns       Kruskal-Wallis
## 7 Verrucomicrobia Abundance 0.142  0.5   0.14     ns       Kruskal-Wallis
Means<-compare_means(Abundance ~ Sampling_station, data = dfScrapers, group.by = "Phylum", p.adjust.method = "fdr")

keeps <- c("Phylum","group1","group2","p","method","p.signif")
keeps=Means[keeps]

test3 <- list('Phylum'= keeps$Phylum,'group1'=keeps$group1,'group2'= keeps$group2,'p'=keeps$p)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p>0.1),]            
FilteredResults
##           Phylum group1 group2          p
## 3 Planctomycetes Ostana    PDR 0.06619258
ggplot(dfShredders, aes(x=Sampling_station,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Phylum)+ylab("Relative Abundance Shredders(>1%)")

ggplot(dfScrapers, aes(x=Sampling_station,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Phylum)+ylab("Relative Abundance Scrapers (>1%)")

############
df <- psmelt(FamilyLevel)
df$Abundance=df$Abundance*100
dfShredders<-subset(df,FFG =="Shredder")
dfScrapers<-subset(df,FFG == "Scraper")

Means<-compare_means(Abundance ~ Sampling_station, data = dfShredders, group.by = "Family", p.adjust.method = "fdr")

keeps <- c("Family","group1","group2","p","method","p.signif")
keeps=Means[keeps]

test3 <- list('Family'= keeps$Family,'group1'=keeps$group1,'group2'= keeps$group2,'p'=keeps$p)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p>0.1),]            
FilteredResults
##               Family group1 group2          p
## 3    Lachnospiraceae Ostana    PDR 0.05182993
## 4      Rikenellaceae Ostana    PDR 0.05182993
## 6   Rhodobacteraceae Ostana    PDR 0.05182993
## 7 Enterobacteriaceae Ostana    PDR 0.05182993
print("Scrapers")
## [1] "Scrapers"
compare_means(Abundance ~ Sampling_station, data = dfScrapers, group.by = "Family", p.adjust.method = "fdr",method="kruskal.test")
## # A tibble: 10 x 7
##    Family           .y.             p   p.adj p.format p.signif method     
##    <fct>            <chr>       <dbl>   <dbl> <chr>    <chr>    <chr>      
##  1 Enterobacteriac~ Abundan~   0.0143   0.086 0.014    *        Kruskal-Wa~
##  2 Aeromonadaceae   Abundan~   0.283    0.37  0.283    ns       Kruskal-Wa~
##  3 Rhodobacteraceae Abundan~   0.142    0.34  0.142    ns       Kruskal-Wa~
##  4 Comamonadaceae   Abundan~   0.169    0.34  0.169    ns       Kruskal-Wa~
##  5 Methylophilaceae Abundan~   0.371    0.37  0.371    ns       Kruskal-Wa~
##  6 Cytophagaceae    Abundan~   0.371    0.37  0.371    ns       Kruskal-Wa~
##  7 Desulfovibriona~ Abundan~ NaN      NaN     NA       ?        Kruskal-Wa~
##  8 Rikenellaceae    Abundan~ NaN      NaN     NA       ?        Kruskal-Wa~
##  9 Ruminococcaceae  Abundan~ NaN      NaN     NA       ?        Kruskal-Wa~
## 10 Lachnospiraceae  Abundan~ NaN      NaN     NA       ?        Kruskal-Wa~
Means<-compare_means(Abundance ~ Sampling_station, data = dfScrapers, group.by = "Family", p.adjust.method = "fdr")

keeps <- c("Family","group1","group2","p","method","p.signif")
keeps=Means[keeps]

test3 <- list('Family'= keeps$Family,'group1'=keeps$group1,'group2'= keeps$group2,'p'=keeps$p)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p>0.1),]            
FilteredResults
##               Family group1 group2          p
## 1 Enterobacteriaceae Ostana    PDR 0.01996445
ggplot(dfShredders, aes(x=Sampling_station,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Family)+ylab("Relative Abundance Shredders(>1%)")

ggplot(dfScrapers, aes(x=Sampling_station,y=Abundance,fill=FFG))+geom_boxplot()+ scale_fill_manual(values=cbPalette)+facet_wrap(~Family)+ylab("Relative Abundance Scrapers (>1%)")

Family Level

The top 5 Families across all samples are given in the table below

df <- psmelt(FamilyLevel)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Family", "FFG"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
dfStats<-subset(df,FFG!="Predator")
dfStats<-subset(dfStats,FFG!="Filterer")
df2 <- psmelt(FamilyAll)
df2$Abundance<-df2$Abundance*100
FamilyAll<-ddply(df2, c("Family"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
FamilySorted<-FamilyAll[order(-FamilyAll$mean),]
#FamilySorted[1:5,]

FamilyAll2<-ddply(df2, c("Family","FFG"), summarise, #For Rel abu tabs
                  N    = length(Abundance),
                  mean = mean(Abundance),
                  sd   = sd(Abundance),
                  se   = sd / sqrt(N)
)
FamilySorted2<-FamilyAll2[order(-FamilyAll$mean),]
FamilySorted2
##                    Family      FFG N        mean          sd          se
## 62         Aeromonadaceae Predator 6 0.000000000  0.00000000 0.000000000
## 56      Acidobacteriaceae Shredder 7 0.004761905  0.01259882 0.004761905
## 126     Brevibacteriaceae Predator 6 0.038888889  0.09525793 0.038888889
## 43                  A1-B1  Scraper 9 0.000000000  0.00000000 0.000000000
## 49       Acetobacteraceae Filterer 4 0.050000000  0.05773503 0.028867513
## 133      Burkholderiaceae Filterer 4 0.000000000  0.00000000 0.000000000
## 84     Anaeroplasmataceae Shredder 7 0.090476190  0.08758971 0.033105799
## 16    [Exiguobacteraceae] Shredder 7 0.000000000  0.00000000 0.000000000
## 132          Brucellaceae Shredder 7 0.014285714  0.02622653 0.009912695
## 90       Armatimonadaceae Predator 6 0.155555556  0.15587269 0.063634760
## 135      Burkholderiaceae  Scraper 9 0.029629630  0.08888889 0.029629630
## 119    Bifidobacteriaceae  Scraper 9 0.003703704  0.01111111 0.003703704
## 149     Chamaesiphonaceae Filterer 4 0.000000000  0.00000000 0.000000000
## 95              auto67_4W  Scraper 9 0.000000000  0.00000000 0.000000000
## 33              0319-6G20 Filterer 4 0.000000000  0.00000000 0.000000000
## 39                211ds20  Scraper 9 0.000000000  0.00000000 0.000000000
## 69         Alcaligenaceae Filterer 4 0.000000000  0.00000000 0.000000000
## 147      Caulobacteraceae  Scraper 9 0.125925926  0.37777778 0.125925926
## 89       Armatimonadaceae Filterer 4 0.016666667  0.01924501 0.009622504
## 117    Bifidobacteriaceae Filterer 4 0.000000000  0.00000000 0.000000000
## 50       Acetobacteraceae Predator 6 0.016666667  0.02788867 0.011385501
## 114      Beijerinckiaceae Predator 6 0.000000000  0.00000000 0.000000000
## 150     Chamaesiphonaceae Predator 6 0.027777778  0.05340273 0.021801574
## 141 Caldicoprobacteraceae Filterer 4 0.000000000  0.00000000 0.000000000
## 98            Bacillaceae Predator 6 0.011111111  0.02721655 0.011111111
## 78       Anaerobrancaceae Predator 6 0.000000000  0.00000000 0.000000000
## 115      Beijerinckiaceae  Scraper 9 0.400000000  0.80484643 0.268282144
## 146      Caulobacteraceae Predator 6 0.277777778  0.13109228 0.053518198
## 113      Beijerinckiaceae Filterer 4 0.000000000  0.00000000 0.000000000
## 128     Brevibacteriaceae Shredder 7 0.000000000  0.00000000 0.000000000
## 54      Acidobacteriaceae Predator 6 0.000000000  0.00000000 0.000000000
## 123     Bradyrhizobiaceae  Scraper 9 0.122222222  0.25927249 0.086424162
## 73       Alteromonadaceae Filterer 4 0.000000000  0.00000000 0.000000000
## 129          Brucellaceae Filterer 4 0.008333333  0.01666667 0.008333333
## 86        Anaplasmataceae Predator 6 0.027777778  0.06804138 0.027777778
## 3       [Bryobacteraceae]  Scraper 9 0.000000000  0.00000000 0.000000000
## 104    Bacteriovoracaceae Shredder 7 0.000000000  0.00000000 0.000000000
## 64         Aeromonadaceae Shredder 7 0.000000000  0.00000000 0.000000000
## 29        [Weeksellaceae] Filterer 4 0.008333333  0.01666667 0.008333333
## 48                    A4b Shredder 7 0.000000000  0.00000000 0.000000000
## 37                211ds20 Filterer 4 0.000000000  0.00000000 0.000000000
## 41                  A1-B1 Filterer 4 0.000000000  0.00000000 0.000000000
## 28      [Tissierellaceae] Shredder 7 0.000000000  0.00000000 0.000000000
## 42                  A1-B1 Predator 6 0.038888889  0.06116160 0.024969117
## 106        Bacteroidaceae Predator 6 0.000000000  0.00000000 0.000000000
## 125     Brevibacteriaceae Filterer 4 0.000000000  0.00000000 0.000000000
## 61         Aeromonadaceae Filterer 4 0.000000000  0.00000000 0.000000000
## 75       Alteromonadaceae  Scraper 9 0.000000000  0.00000000 0.000000000
## 6       [Cerasicoccaceae] Predator 6 0.005555556  0.01360828 0.005555556
## 8       [Cerasicoccaceae] Shredder 7 0.000000000  0.00000000 0.000000000
## 140                  C111 Shredder 7 0.019047619  0.01781742 0.006734350
## 25      [Tissierellaceae] Filterer 4 0.000000000  0.00000000 0.000000000
## 107        Bacteroidaceae  Scraper 9 0.000000000  0.00000000 0.000000000
## 143 Caldicoprobacteraceae  Scraper 9 0.000000000  0.00000000 0.000000000
## 118    Bifidobacteriaceae Predator 6 0.000000000  0.00000000 0.000000000
## 57       Actinomycetaceae Filterer 4 0.000000000  0.00000000 0.000000000
## 124     Bradyrhizobiaceae Shredder 7 0.000000000  0.00000000 0.000000000
## 138                  C111 Predator 6 0.150000000  0.16295875 0.066527633
## 68             AK1AB1_02E Shredder 7 0.000000000  0.00000000 0.000000000
## 92       Armatimonadaceae Shredder 7 0.095238095  0.11126973 0.042056004
## 26      [Tissierellaceae] Predator 6 0.000000000  0.00000000 0.000000000
## 35              0319-6G20  Scraper 9 0.000000000  0.00000000 0.000000000
## 23     [Mogibacteriaceae]  Scraper 9 0.000000000  0.00000000 0.000000000
## 102    Bacteriovoracaceae Predator 6 0.283333333  0.14414499 0.058846945
## 32        [Weeksellaceae] Shredder 7 0.085714286  0.09786072 0.036987874
## 45                    A4b Filterer 4 0.000000000  0.00000000 0.000000000
## 139                  C111  Scraper 9 0.040740741  0.12222222 0.040740741
## 13    [Exiguobacteraceae] Filterer 4 0.000000000  0.00000000 0.000000000
## 20    [Fimbriimonadaceae] Shredder 7 0.009523810  0.01626500 0.006147593
## 31        [Weeksellaceae]  Scraper 9 0.000000000  0.00000000 0.000000000
## 24     [Mogibacteriaceae] Shredder 7 0.352380952  0.36201961 0.136830553
## 44                  A1-B1 Shredder 7 0.042857143  0.07382232 0.027902216
## 58       Actinomycetaceae Predator 6 0.000000000  0.00000000 0.000000000
## 101    Bacteriovoracaceae Filterer 4 0.016666667  0.03333333 0.016666667
## 105        Bacteroidaceae Filterer 4 0.000000000  0.00000000 0.000000000
## 144 Caldicoprobacteraceae Shredder 7 0.028571429  0.04049953 0.015307382
## 19    [Fimbriimonadaceae]  Scraper 9 0.000000000  0.00000000 0.000000000
## 91       Armatimonadaceae  Scraper 9 0.000000000  0.00000000 0.000000000
## 110    Bdellovibrionaceae Predator 6 0.244444444  0.22377237 0.091354688
## 122     Bradyrhizobiaceae Predator 6 0.005555556  0.01360828 0.005555556
## 53      Acidobacteriaceae Filterer 4 0.000000000  0.00000000 0.000000000
## 38                211ds20 Predator 6 0.005555556  0.01360828 0.005555556
## 2       [Bryobacteraceae] Predator 6 0.083333333  0.11105554 0.045338235
## 5       [Cerasicoccaceae] Filterer 4 0.000000000  0.00000000 0.000000000
## 134      Burkholderiaceae Predator 6 0.000000000  0.00000000 0.000000000
## 77       Anaerobrancaceae Filterer 4 0.000000000  0.00000000 0.000000000
## 21     [Mogibacteriaceae] Filterer 4 0.025000000  0.03191424 0.015957118
## 52       Acetobacteraceae Shredder 7 0.023809524  0.04178554 0.015793451
## 103    Bacteriovoracaceae  Scraper 9 0.000000000  0.00000000 0.000000000
## 99            Bacillaceae  Scraper 9 0.237037037  0.71111111 0.237037037
## 81     Anaeroplasmataceae Filterer 4 0.000000000  0.00000000 0.000000000
## 142 Caldicoprobacteraceae Predator 6 0.000000000  0.00000000 0.000000000
## 11  [Chthoniobacteraceae]  Scraper 9 0.166666667  0.28037673 0.093458910
## 130          Brucellaceae Predator 6 0.066666667  0.06666667 0.027216553
## 1       [Bryobacteraceae] Filterer 4 0.000000000  0.00000000 0.000000000
## 72         Alcaligenaceae Shredder 7 0.000000000  0.00000000 0.000000000
## 120    Bifidobacteriaceae Shredder 7 0.000000000  0.00000000 0.000000000
## 46                    A4b Predator 6 0.022222222  0.04036867 0.016480441
## 63         Aeromonadaceae  Scraper 9 8.733333333 15.74738955 5.249129851
## 74       Alteromonadaceae Predator 6 0.111111111  0.22673202 0.092562958
## 97            Bacillaceae Filterer 4 0.000000000  0.00000000 0.000000000
## 14    [Exiguobacteraceae] Predator 6 0.005555556  0.01360828 0.005555556
## 85        Anaplasmataceae Filterer 4 0.000000000  0.00000000 0.000000000
## 112    Bdellovibrionaceae Shredder 7 0.028571429  0.05245305 0.019825390
## 79       Anaerobrancaceae  Scraper 9 0.000000000  0.00000000 0.000000000
## 100           Bacillaceae Shredder 7 0.009523810  0.02519763 0.009523810
## 18    [Fimbriimonadaceae] Predator 6 0.100000000  0.11737878 0.047919686
## 34              0319-6G20 Predator 6 0.016666667  0.04082483 0.016666667
## 65             AK1AB1_02E Filterer 4 0.008333333  0.01666667 0.008333333
## 36              0319-6G20 Shredder 7 0.000000000  0.00000000 0.000000000
## 121     Bradyrhizobiaceae Filterer 4 0.000000000  0.00000000 0.000000000
## 22     [Mogibacteriaceae] Predator 6 0.000000000  0.00000000 0.000000000
## 83     Anaeroplasmataceae  Scraper 9 0.000000000  0.00000000 0.000000000
## 82     Anaeroplasmataceae Predator 6 0.000000000  0.00000000 0.000000000
## 94              auto67_4W Predator 6 0.150000000  0.16158933 0.065968567
## 131          Brucellaceae  Scraper 9 4.648148148  6.95218014 2.317393379
## 12  [Chthoniobacteraceae] Shredder 7 0.038095238  0.06214848 0.023489918
## 59       Actinomycetaceae  Scraper 9 0.000000000  0.00000000 0.000000000
## 109    Bdellovibrionaceae Filterer 4 0.158333333  0.14240006 0.071200031
## 9   [Chthoniobacteraceae] Filterer 4 0.275000000  0.18534253 0.092671263
## 66             AK1AB1_02E Predator 6 0.000000000  0.00000000 0.000000000
## 67             AK1AB1_02E  Scraper 9 0.000000000  0.00000000 0.000000000
## 127     Brevibacteriaceae  Scraper 9 0.129629630  0.38888889 0.129629630
## 148      Caulobacteraceae Shredder 7 0.085714286  0.22677868 0.085714286
## 88        Anaplasmataceae Shredder 7 0.000000000  0.00000000 0.000000000
## 145      Caulobacteraceae Filterer 4 0.016666667  0.03333333 0.016666667
## 15    [Exiguobacteraceae]  Scraper 9 0.000000000  0.00000000 0.000000000
## 40                211ds20 Shredder 7 0.000000000  0.00000000 0.000000000
## 55      Acidobacteriaceae  Scraper 9 0.037037037  0.11111111 0.037037037
## 60       Actinomycetaceae Shredder 7 0.009523810  0.02519763 0.009523810
## 87        Anaplasmataceae  Scraper 9 0.000000000  0.00000000 0.000000000
## 93              auto67_4W Filterer 4 0.050000000  0.05773503 0.028867513
## 96              auto67_4W Shredder 7 0.000000000  0.00000000 0.000000000
## 108        Bacteroidaceae Shredder 7 0.004761905  0.01259882 0.004761905
## 111    Bdellovibrionaceae  Scraper 9 0.114814815  0.23339946 0.077799821
## 136      Burkholderiaceae Shredder 7 0.000000000  0.00000000 0.000000000
## 137                  C111 Filterer 4 0.091666667  0.18333333 0.091666667
## 4       [Bryobacteraceae] Shredder 7 0.000000000  0.00000000 0.000000000
## 7       [Cerasicoccaceae]  Scraper 9 0.081481481  0.18189673 0.060632243
## 10  [Chthoniobacteraceae] Predator 6 0.333333333  0.36757463 0.150061716
## 17    [Fimbriimonadaceae] Filterer 4 0.016666667  0.03333333 0.016666667
## 27      [Tissierellaceae]  Scraper 9 0.003703704  0.01111111 0.003703704
## 30        [Weeksellaceae] Predator 6 0.300000000  0.35402448 0.144529889
## 47                    A4b  Scraper 9 0.000000000  0.00000000 0.000000000
## 51       Acetobacteraceae  Scraper 9 0.088888889  0.16914819 0.056382731
## 70         Alcaligenaceae Predator 6 0.000000000  0.00000000 0.000000000
## 71         Alcaligenaceae  Scraper 9 0.029629630  0.07718024 0.025726748
## 76       Alteromonadaceae Shredder 7 0.042857143  0.11338934 0.042857143
## 80       Anaerobrancaceae Shredder 7 0.161904762  0.30088054 0.113722155
## 116      Beijerinckiaceae Shredder 7 0.000000000  0.00000000 0.000000000
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = Family),colour="black", stat="identity")+xlab("FFG")+ylab("Relative Abundance (> 2%)")+theme(axis.title.x=element_blank())+ theme(axis.text.x = element_text(angle = 45, hjust = 1))+ scale_fill_manual(values=cbPalette)
cdataplot#+  scale_x_discrete(labels=c("0 hrs", "24 hrs", "48 hrs","72 hrs"))

#position = "fill"#+geom_errorbar(aes(ymin=mean-se,ymax=mean+se),color="black",width=1)+geom_line(size=1.5,linetype="dashed")+geom_point(size=6)+ylab(

Trtdata2 <- ddply(df, c("Family", "FFG","Sampling_station"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
cdataplot2=ggplot(Trtdata2, aes(x=FFG,y=mean))+geom_bar(aes(fill = Family),colour="black", stat="identity")+xlab("FFG")+ylab("Relative Abundance (> 1%)")+theme(axis.title.x=element_blank())+ theme(axis.text.x = element_text(angle = 45, hjust = 1))+ scale_fill_manual(values=cbPalette)+facet_wrap(~Sampling_station)
cdataplot2

########
dfStats<-subset(df, FFG!="Predator")
dfStats<-subset(dfStats,FFG!="Filterer")
compare_means(Abundance ~ FFG, data = dfStats, group.by = "Family", p.adjust.method = "fdr",method="kruskal.test")# Differences between shredder and scrapers only
## # A tibble: 10 x 7
##    Family           .y.             p   p.adj p.format p.signif method     
##    <fct>            <chr>       <dbl>   <dbl> <chr>    <chr>    <chr>      
##  1 Enterobacteriac~ Abundan~ 0.340    3.80e-1 0.34041  ns       Kruskal-Wa~
##  2 Aeromonadaceae   Abundan~ 0.0516   6.40e-2 0.05155  ns       Kruskal-Wa~
##  3 Rhodobacteraceae Abundan~ 0.958    9.60e-1 0.95779  ns       Kruskal-Wa~
##  4 Ruminococcaceae  Abundan~ 0.000239 6.00e-4 0.00024  ***      Kruskal-Wa~
##  5 Desulfovibriona~ Abundan~ 0.000239 6.00e-4 0.00024  ***      Kruskal-Wa~
##  6 Lachnospiraceae  Abundan~ 0.000239 6.00e-4 0.00024  ***      Kruskal-Wa~
##  7 Rikenellaceae    Abundan~ 0.000239 6.00e-4 0.00024  ***      Kruskal-Wa~
##  8 Methylophilaceae Abundan~ 0.000369 6.10e-4 0.00037  ***      Kruskal-Wa~
##  9 Comamonadaceae   Abundan~ 0.000818 1.20e-3 0.00082  ***      Kruskal-Wa~
## 10 Cytophagaceae    Abundan~ 0.000369 6.10e-4 0.00037  ***      Kruskal-Wa~
#df<-subset(df, Family!="Aeromonadaceae"&Family!="Rhodobacteraceae"&Family!="Enterobacteriaceae")# Remove groups that were not sig by Kruksal-Wallis

#df <- df[which(df$Family!="Aeromonadaceae")] #& df$Family!="Rhodobacteraceae"& df$Family!="Enterobacteriaceae"
#df
Means<-compare_means(Abundance ~ FFG, data = df, group.by = "Family", p.adjust.method = "fdr")


Trtdata <- ddply(df, c("Family", "FFG"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
SigList<-length(unique(Trtdata$Family))
for (i in levels(Means$Family)){
  Tax<-i
  TaxAbundance<-subset(Means,Family==i )
  Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
  difference<-TaxAbundance$p
  names(difference)<-Hyphenated
  Letters<-multcompLetters(difference)
  #print(Letters)
  SigList[i]<-Letters
  
}
vec<-unlist(SigList)
vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Family)+xlab("FFG")+ylab("Relative Abundance (> 2%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Family)+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+geom_text(aes(x=FFG, y=mean+se+5,label=vec))+scale_fill_manual(values=cbPalette)
cdataplot

(Trtdata)
##                 Family      FFG N         mean          sd           se
## 1       Aeromonadaceae Filterer 4  0.000000000  0.00000000  0.000000000
## 2       Aeromonadaceae Predator 6  0.000000000  0.00000000  0.000000000
## 3       Aeromonadaceae  Scraper 9  8.733333333 15.74738955  5.249129851
## 4       Aeromonadaceae Shredder 7  0.000000000  0.00000000  0.000000000
## 5       Comamonadaceae Filterer 4 12.333333333  8.80151502  4.400757511
## 6       Comamonadaceae Predator 6  6.294444444  2.43715833  0.994965723
## 7       Comamonadaceae  Scraper 9  0.362962963  0.46050831  0.153502769
## 8       Comamonadaceae Shredder 7  2.314285714  1.02068553  0.385782867
## 9        Cytophagaceae Filterer 4  1.500000000  1.09273696  0.546368482
## 10       Cytophagaceae Predator 6 14.833333333  8.29572849  3.386716973
## 11       Cytophagaceae  Scraper 9  0.011111111  0.03333333  0.011111111
## 12       Cytophagaceae Shredder 7  1.171428571  0.61626971  0.232928057
## 13 Desulfovibrionaceae Filterer 4 10.708333333  5.08821259  2.544106297
## 14 Desulfovibrionaceae Predator 6  0.255555556  0.62598071  0.255555556
## 15 Desulfovibrionaceae  Scraper 9  0.000000000  0.00000000  0.000000000
## 16 Desulfovibrionaceae Shredder 7  9.923809524  3.25671468  1.230922446
## 17  Enterobacteriaceae Filterer 4  0.025000000  0.05000000  0.025000000
## 18  Enterobacteriaceae Predator 6  8.122222222 15.38135475  6.279411783
## 19  Enterobacteriaceae  Scraper 9 24.911111111 39.25302182 13.084340608
## 20  Enterobacteriaceae Shredder 7  1.985714286  2.21926390  0.838802912
## 21     Lachnospiraceae Filterer 4  6.000000000  5.05253878  2.526269391
## 22     Lachnospiraceae Predator 6  0.011111111  0.02721655  0.011111111
## 23     Lachnospiraceae  Scraper 9  0.000000000  0.00000000  0.000000000
## 24     Lachnospiraceae Shredder 7  7.909523810  3.33042730  1.258783201
## 25    Methylophilaceae Filterer 4  2.508333333  1.45178689  0.725893447
## 26    Methylophilaceae Predator 6  3.155555556  3.02571693  1.235243766
## 27    Methylophilaceae  Scraper 9  0.011111111  0.03333333  0.011111111
## 28    Methylophilaceae Shredder 7  5.090476190  2.19145768  0.828293148
## 29    Rhodobacteraceae Filterer 4  3.750000000  3.46297881  1.731489404
## 30    Rhodobacteraceae Predator 6  7.327777778  3.02591890  1.235326218
## 31    Rhodobacteraceae  Scraper 9  3.977777778  6.42209727  2.140699090
## 32    Rhodobacteraceae Shredder 7  2.647619048  2.44198289  0.922982775
## 33       Rikenellaceae Filterer 4  4.725000000  4.26287809  2.131439046
## 34       Rikenellaceae Predator 6  0.000000000  0.00000000  0.000000000
## 35       Rikenellaceae  Scraper 9  0.000000000  0.00000000  0.000000000
## 36       Rikenellaceae Shredder 7  7.076190476  3.38677790  1.280081725
## 37     Ruminococcaceae Filterer 4  5.250000000  3.50518135  1.752590675
## 38     Ruminococcaceae Predator 6  0.005555556  0.01360828  0.005555556
## 39     Ruminococcaceae  Scraper 9  0.000000000  0.00000000  0.000000000
## 40     Ruminococcaceae Shredder 7  9.690476190  4.86722654  1.839638714
TrtdataPlot<-subset(Trtdata, Family!="Aeromonadaceae"&Family!="Rhodobacteraceae"&Family!="Enterobacteriaceae"&Family!="Cytophagaceae")
(TrtdataPlot)
##                 Family      FFG N         mean         sd          se
## 5       Comamonadaceae Filterer 4 12.333333333 8.80151502 4.400757511
## 6       Comamonadaceae Predator 6  6.294444444 2.43715833 0.994965723
## 7       Comamonadaceae  Scraper 9  0.362962963 0.46050831 0.153502769
## 8       Comamonadaceae Shredder 7  2.314285714 1.02068553 0.385782867
## 13 Desulfovibrionaceae Filterer 4 10.708333333 5.08821259 2.544106297
## 14 Desulfovibrionaceae Predator 6  0.255555556 0.62598071 0.255555556
## 15 Desulfovibrionaceae  Scraper 9  0.000000000 0.00000000 0.000000000
## 16 Desulfovibrionaceae Shredder 7  9.923809524 3.25671468 1.230922446
## 21     Lachnospiraceae Filterer 4  6.000000000 5.05253878 2.526269391
## 22     Lachnospiraceae Predator 6  0.011111111 0.02721655 0.011111111
## 23     Lachnospiraceae  Scraper 9  0.000000000 0.00000000 0.000000000
## 24     Lachnospiraceae Shredder 7  7.909523810 3.33042730 1.258783201
## 25    Methylophilaceae Filterer 4  2.508333333 1.45178689 0.725893447
## 26    Methylophilaceae Predator 6  3.155555556 3.02571693 1.235243766
## 27    Methylophilaceae  Scraper 9  0.011111111 0.03333333 0.011111111
## 28    Methylophilaceae Shredder 7  5.090476190 2.19145768 0.828293148
## 33       Rikenellaceae Filterer 4  4.725000000 4.26287809 2.131439046
## 34       Rikenellaceae Predator 6  0.000000000 0.00000000 0.000000000
## 35       Rikenellaceae  Scraper 9  0.000000000 0.00000000 0.000000000
## 36       Rikenellaceae Shredder 7  7.076190476 3.38677790 1.280081725
## 37     Ruminococcaceae Filterer 4  5.250000000 3.50518135 1.752590675
## 38     Ruminococcaceae Predator 6  0.005555556 0.01360828 0.005555556
## 39     Ruminococcaceae  Scraper 9  0.000000000 0.00000000 0.000000000
## 40     Ruminococcaceae Shredder 7  9.690476190 4.86722654 1.839638714
#vec<-c("a","b","a","b","a","b","a","b","a","b","a","b")
cdataplot=ggplot(TrtdataPlot, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Family)+xlab("FFG")+ylab("Relative Abundance (> 2%)") + theme(axis.title.x=element_blank())+facet_wrap(~Family)+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)#theme(axis.text.x = element_text(angle = 45, hjust = 1))+
cdataplot

MeansStation=compare_means(Abundance ~ Sampling_station, data = dfStats, group.by = "Family", p.adjust.method = "fdr")


#NComparisons<-length(unique(metadata$FFG))*length(unique(Trtdata$Family))
#SigList<-length(unique(Trtdata$Family))
#SigLetters2<-vector(length=NComparisons)
#vec<-unlist(lst)
#Means=compare_means(Abundance ~ FFG, data = dfStats, group.by = "Family", p.adjust.method = "fdr")
#for (i in levels(Means$Family)){
#  Tax<-i
#  TaxAbundance<-subset(Means,Family==i )
#  Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
#  difference<-TaxAbundance$p.adj
#  names(difference)<-Hyphenated
#  Letters<-multcompLetters(difference)
  #print(Letters)
#  SigList[i]<-Letters
  
#}
#vec<-unlist(SigList)
#vec<-vec[-1]

cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Family)+xlab("FFG")+ylab("Relative Abundance (> 2%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Family)+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
cdataplot

#+geom_text(aes(x=FFG, y=mean+se+5,label=vec))

Means=compare_means(Abundance ~ FFG, data = dfStats, 
                    group.by = "Family", p.adjust.method = "fdr")
keeps <- c("Family","group1","group2","p.adj","method","p.signif")
keeps=Means[keeps]
#keeps
 
 
test3 <- list('Family'= keeps$Family,'group1'=keeps$group1,'group2'= keeps$group2,'p.adj'=keeps$p.adj)
test3= as.data.frame(test3)
FilteredResults<-test3[!(test3$p.adj>0.1),]            
FilteredResults
##                 Family  group1   group2   p.adj
## 2       Aeromonadaceae Scraper Shredder 0.07600
## 4      Ruminococcaceae Scraper Shredder 0.00075
## 5  Desulfovibrionaceae Scraper Shredder 0.00075
## 6      Lachnospiraceae Scraper Shredder 0.00075
## 7        Rikenellaceae Scraper Shredder 0.00075
## 8     Methylophilaceae Scraper Shredder 0.00076
## 9       Comamonadaceae Scraper Shredder 0.00140
## 10       Cytophagaceae Scraper Shredder 0.00076

Genus level

The top 5 genera across all samples are given in the table below

df <- psmelt(GenusLevel)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Genus", "FFG"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
df2 <- psmelt(GenusAll)
df2$Abundance<-df2$Abundance*100
GenusAll<-ddply(df2, c("Genus"), summarise,
                N    = length(Abundance),
                mean = mean(Abundance),
                sd   = sd(Abundance),
                se   = sd / sqrt(N)
)
GenusSorted<-GenusAll[order(-GenusAll$mean),]
GenusSorted[1:5,]
##              Genus  N     mean       sd        se
## 131    Rhodobacter 26 4.178205 4.430901 0.8689712
## 101  Methylotenera 26 2.488462 2.729264 0.5352528
## 88  Leadbetterella 26 1.992308 5.600399 1.0983287
## 58    Dysgonomonas 26 1.851282 3.369883 0.6608884
## 40     Clostridium 52 1.410897 3.143557 0.4359329
GenusAll2<-ddply(df2, c("Genus","FFG"), summarise, #For Rel abu tabs
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
GenusSorted2<-GenusAll2[order(-GenusAll$mean),]
dfStats<-subset(df,FFG!="Predator")
dfStats<-subset(dfStats,FFG!="Filterer")
Trtdata2 <- ddply(df, c("Genus", "FFG","Sampling_station"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
cdataplot2=ggplot(Trtdata2, aes(x=FFG,y=mean))+geom_bar(aes(fill = Genus),colour="black", stat="identity")+xlab("FFG")+ylab("Relative Abundance (> 1%)")+theme(axis.title.x=element_blank())+ theme(axis.text.x = element_text(angle = 45, hjust = 1))+ scale_fill_manual(values=cbPalette)+facet_wrap(~Sampling_station)
cdataplot2

#NComparisons<-length(unique(metadata$FFG))*length(unique(Trtdata$Genus))
#SigList<-length(unique(Trtdata$Genus))
#SigLetters2<-vector(length=NComparisons)

Means=compare_means(Abundance ~ FFG, data = dfStats, group.by = "Genus", p.adjust.method = "fdr")
MeansStation=compare_means(Abundance ~ Sampling_station, data = dfStats, group.by = "Genus", p.adjust.method = "fdr")
#vec<-unlist(lst)
#for (i in levels(Means$Genus)){
#  Tax<-i
#  TaxAbundance<-subset(Means,Genus==i )
  #print(TaxAbundance)
#  Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
#  difference<-TaxAbundance$p.adj
#  names(difference)<-Hyphenated
#  Letters<-multcompLetters(difference)
  #print(Letters)
#  SigList[i]<-Letters
  
#}
#vec<-unlist(SigList)
#vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = Genus),colour="black", stat="identity")+xlab("FFG")+ylab("Relative Abundance (> 2%)") + theme(axis.text.x = element_text(angle = 0, hjust = 0.5))+theme(axis.title.x=element_blank())+ theme(axis.text.x = element_text(angle = 45, hjust = 1))+ scale_fill_manual(values=cbPalette)#+geom_text(aes(x=FFG, y=mean+se+10,label=vec))
cdataplot#+  scale_x_discrete(labels=c("0 hrs", "24 hrs", "48 hrs","72 hrs"))

cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Genus)+xlab("FFG")+ylab("Relative Abundance (> 2%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Genus)+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
#position = "fill"#+geom_errorbar(aes(ymin=mean-se,ymax=mean+se),color="black",width=1)+geom_line(size=1.5,linetype="dashed")+geom_point(size=6)+ylab(geom_text(aes(x=FFG, y=mean+se+5,label=vec))

#+ theme(axis.text.x = element_text(angle = 45, hjust = 1))#+ scale_fill_manual(values=cbPalette)
compare_means(Abundance~FFG, data=dfStats,group.by="Genus",method = "kruskal.test",p.adjust.method="fdr")
## # A tibble: 3 x 7
##   Genus         .y.              p   p.adj p.format p.signif method        
##   <fct>         <chr>        <dbl>   <dbl> <chr>    <chr>    <chr>         
## 1 Rhodobacter   Abundance 0.958    0.96    0.95776  ns       Kruskal-Wallis
## 2 Clostridium   Abundance 0.000239 0.00055 0.00024  ***      Kruskal-Wallis
## 3 Methylotenera Abundance 0.000369 0.00055 0.00037  ***      Kruskal-Wallis
Means=compare_means(Abundance ~ FFG, data = dfStats, 
                    group.by = "Genus", p.adjust.method = "fdr")
Means
## # A tibble: 3 x 9
##   Genus     .y.     group1  group2       p   p.adj p.format p.signif method
##   <fct>     <chr>   <chr>   <chr>    <dbl>   <dbl> <chr>    <chr>    <chr> 
## 1 Rhodobac~ Abunda~ Scraper Shred~ 1.00e+0 1       1.00000  ns       Wilco~
## 2 Clostrid~ Abunda~ Scraper Shred~ 2.99e-4 0.00068 0.00030  ***      Wilco~
## 3 Methylot~ Abunda~ Scraper Shred~ 4.57e-4 0.00068 0.00046  ***      Wilco~
MeansStation=compare_means(Abundance ~ Sampling_station, data = dfStats, 
                    group.by = "Genus", p.adjust.method = "fdr")
MeansStation
## # A tibble: 3 x 9
##   Genus       .y.      group1 group2      p p.adj p.format p.signif method 
##   <fct>       <chr>    <chr>  <chr>   <dbl> <dbl> <chr>    <chr>    <chr>  
## 1 Rhodobacter Abundan~ Ostana PDR    0.0135 0.041 0.014    *        Wilcox~
## 2 Clostridium Abundan~ Ostana PDR    0.325  0.49  0.325    ns       Wilcox~
## 3 Methyloten~ Abundan~ Ostana PDR    0.866  0.87  0.866    ns       Wilcox~
# #head(Means)
# keeps <- c("Genus","group1","group2","p.adj","method","p.signif")
# keeps=Means[keeps]
# #keeps
# 
# 
# test3 <- list('Genus'= keeps$Genus,'group1'=keeps$group1,'group2'= keeps$group2,'p.adj'=keeps$p.adj)
# test3= as.data.frame(test3)
# #test3
# FilteredResults<-test3[!(test3$p.adj>0.05),]            
# FilteredResults

PhylumAbuAll

Below is a list of the relative abundance of all taxa at the Phylum level by FFG

df <- psmelt(PhylumLevel)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Phylum", "FFG"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
compare_means(Abundance~FFG, data=df,group.by="Phylum",method = "kruskal.test",p.adjust.method="fdr")
## # A tibble: 7 x 7
##   Phylum          .y.              p  p.adj p.format p.signif method       
##   <fct>           <chr>        <dbl>  <dbl> <chr>    <chr>    <chr>        
## 1 Proteobacteria  Abundance 0.0147   0.021  0.01466  *        Kruskal-Wall~
## 2 Tenericutes     Abundance 0.0643   0.075  0.06432  ns       Kruskal-Wall~
## 3 Bacteroidetes   Abundance 0.000400 0.0014 0.00040  ***      Kruskal-Wall~
## 4 Firmicutes      Abundance 0.000209 0.0014 0.00021  ***      Kruskal-Wall~
## 5 Planctomycetes  Abundance 0.0127   0.021  0.01274  *        Kruskal-Wall~
## 6 Verrucomicrobia Abundance 0.00184  0.0043 0.00184  **       Kruskal-Wall~
## 7 Actinobacteria  Abundance 0.118    0.12   0.11762  ns       Kruskal-Wall~
Means=compare_means(Abundance ~ FFG, data = df, 
                    group.by = "Phylum", p.adjust.method = "fdr")

Multiple comparisons for all taxa grouped by Phylum

NComparisons<-length(unique(metadata$FFG))*length(unique(Trtdata$Phylum))
SigList<-length(unique(Trtdata$Phylum))
SigLetters2<-vector(length=NComparisons)
#vec<-unlist(lst)
#for (i in levels(Means$Phylum)){
#  Tax<-i
#  TaxAbundance<-subset(Means,Phylum==i )
#  print(TaxAbundance)
#  Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
#  difference<-TaxAbundance$p.adj
#  names(difference)<-Hyphenated
#  Letters<-multcompLetters(difference)
#  print(Letters)
#  
#}

FamilyAbuAll

Below is a list of the relative abundance of all taxa at the Family level by FFG

df <- psmelt(FamilyLevel)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Family", "FFG"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
compare_means(Abundance~FFG, data=df,group.by="Family",method = "kruskal.test",p.adjust.method="fdr")
## # A tibble: 10 x 7
##    Family          .y.              p   p.adj p.format p.signif method     
##    <fct>           <chr>        <dbl>   <dbl> <chr>    <chr>    <chr>      
##  1 Enterobacteria~ Abundan~ 0.0518    5.20e-2 0.05176  ns       Kruskal-Wa~
##  2 Aeromonadaceae  Abundan~ 0.0365    4.60e-2 0.03654  *        Kruskal-Wa~
##  3 Cytophagaceae   Abundan~ 0.0000527 1.50e-4 5.3e-05  ****     Kruskal-Wa~
##  4 Comamonadaceae  Abundan~ 0.0000611 1.50e-4 6.1e-05  ****     Kruskal-Wa~
##  5 Rhodobacterace~ Abundan~ 0.0464    5.20e-2 0.04640  *        Kruskal-Wa~
##  6 Ruminococcaceae Abundan~ 0.0000556 1.50e-4 5.6e-05  ****     Kruskal-Wa~
##  7 Desulfovibrion~ Abundan~ 0.0000686 1.50e-4 6.9e-05  ****     Kruskal-Wa~
##  8 Lachnospiraceae Abundan~ 0.0000763 1.50e-4 7.6e-05  ****     Kruskal-Wa~
##  9 Rikenellaceae   Abundan~ 0.000129  2.10e-4 0.00013  ***      Kruskal-Wa~
## 10 Methylophilace~ Abundan~ 0.000339  4.80e-4 0.00034  ***      Kruskal-Wa~
Means=compare_means(Abundance ~ FFG, data = df, 
                    group.by = "Family", p.adjust.method = "fdr")

Multiple comparisons for all taxa grouped by Family

NComparisons<-length(unique(metadata$FFG))*length(unique(Trtdata$Family))
SigList<-length(unique(Trtdata$Family))
SigLetters2<-vector(length=NComparisons)
#vec<-unlist(lst)
for (i in levels(Means$Family)){
  Tax<-i
  TaxAbundance<-subset(Means,Family==i )
  print(TaxAbundance)
  Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
  difference<-TaxAbundance$p.adj
  names(difference)<-Hyphenated
  Letters<-multcompLetters(difference)
  print(Letters)
  
}
## # A tibble: 3 x 9
##   Family     .y.     group1  group2       p p.adj p.format p.signif method 
##   <fct>      <chr>   <chr>   <chr>    <dbl> <dbl> <chr>    <chr>    <chr>  
## 1 Aeromonad~ Abunda~ Filter~ Scraper 0.158  0.22  0.15751  ns       Wilcox~
## 2 Aeromonad~ Abunda~ Predat~ Scraper 0.0820 0.13  0.08197  ns       Wilcox~
## 3 Aeromonad~ Abunda~ Scraper Shredd~ 0.0605 0.097 0.06048  ns       Wilcox~
## $Letters
## Filterer Predator  Scraper Shredder 
##      "a"      "a"      "a"      "a" 
## 
## $LetterMatrix
##             a
## Filterer TRUE
## Predator TRUE
## Scraper  TRUE
## Shredder TRUE
## 
## # A tibble: 6 x 9
##   Family    .y.     group1  group2        p  p.adj p.format p.signif method
##   <fct>     <chr>   <chr>   <chr>     <dbl>  <dbl> <chr>    <chr>    <chr> 
## 1 Comamona~ Abunda~ Filter~ Predat~ 2.40e-1 0.31   0.23953  ns       Wilco~
## 2 Comamona~ Abunda~ Filter~ Scraper 6.55e-3 0.016  0.00655  **       Wilco~
## 3 Comamona~ Abunda~ Filter~ Shredd~ 1.07e-2 0.023  0.01073  *        Wilco~
## 4 Comamona~ Abunda~ Predat~ Scraper 1.69e-3 0.0073 0.00169  **       Wilco~
## 5 Comamona~ Abunda~ Predat~ Shredd~ 5.28e-3 0.014  0.00528  **       Wilco~
## 6 Comamona~ Abunda~ Scraper Shredd~ 9.89e-4 0.0046 0.00099  ***      Wilco~
## Filterer Predator  Scraper Shredder 
##      "a"      "a"      "b"      "c" 
## # A tibble: 6 x 9
##   Family    .y.     group1  group2        p  p.adj p.format p.signif method
##   <fct>     <chr>   <chr>   <chr>     <dbl>  <dbl> <chr>    <chr>    <chr> 
## 1 Cytophag~ Abunda~ Filter~ Predat~ 1.42e-2 0.028  0.01421  *        Wilco~
## 2 Cytophag~ Abunda~ Filter~ Scraper 2.08e-3 0.0073 0.00208  **       Wilco~
## 3 Cytophag~ Abunda~ Filter~ Shredd~ 9.25e-1 0.94   0.92472  ns       Wilco~
## 4 Cytophag~ Abunda~ Predat~ Scraper 7.06e-4 0.0046 0.00071  ***      Wilco~
## 5 Cytophag~ Abunda~ Predat~ Shredd~ 3.41e-3 0.0095 0.00341  **       Wilco~
## 6 Cytophag~ Abunda~ Scraper Shredd~ 4.57e-4 0.0043 0.00046  ***      Wilco~
## Filterer Predator  Scraper Shredder 
##      "a"      "b"      "c"      "a" 
## # A tibble: 6 x 9
##   Family     .y.     group1  group2       p  p.adj p.format p.signif method
##   <fct>      <chr>   <chr>   <chr>    <dbl>  <dbl> <chr>    <chr>    <chr> 
## 1 Desulfovi~ Abunda~ Filter~ Preda~ 8.91e-3 0.02   0.00891  **       Wilco~
## 2 Desulfovi~ Abunda~ Filter~ Scrap~ 9.73e-4 0.0046 0.00097  ***      Wilco~
## 3 Desulfovi~ Abunda~ Filter~ Shred~ 6.37e-1 0.7    0.63660  ns       Wilco~
## 4 Desulfovi~ Abunda~ Predat~ Scrap~ 2.76e-1 0.34   0.27630  ns       Wilco~
## 5 Desulfovi~ Abunda~ Predat~ Shred~ 2.59e-3 0.0076 0.00259  **       Wilco~
## 6 Desulfovi~ Abunda~ Scraper Shred~ 2.99e-4 0.0042 0.00030  ***      Wilco~
## Filterer Predator  Scraper Shredder 
##      "a"      "b"      "b"      "a" 
## # A tibble: 6 x 9
##   Family      .y.     group1  group2       p p.adj p.format p.signif method
##   <fct>       <chr>   <chr>   <chr>    <dbl> <dbl> <chr>    <chr>    <chr> 
## 1 Enterobact~ Abunda~ Filter~ Predat~ 0.0304 0.052 0.03038  *        Wilco~
## 2 Enterobact~ Abunda~ Filter~ Scraper 0.0189 0.035 0.01892  *        Wilco~
## 3 Enterobact~ Abunda~ Filter~ Shredd~ 0.0334 0.055 0.03336  *        Wilco~
## 4 Enterobact~ Abunda~ Predat~ Scraper 0.768  0.83  0.76828  ns       Wilco~
## 5 Enterobact~ Abunda~ Predat~ Shredd~ 0.886  0.94  0.88625  ns       Wilco~
## 6 Enterobact~ Abunda~ Scraper Shredd~ 0.368  0.44  0.36791  ns       Wilco~
## Filterer Predator  Scraper Shredder 
##      "a"     "ab"      "b"     "ab" 
## # A tibble: 6 x 9
##   Family     .y.     group1  group2       p  p.adj p.format p.signif method
##   <fct>      <chr>   <chr>   <chr>    <dbl>  <dbl> <chr>    <chr>    <chr> 
## 1 Lachnospi~ Abunda~ Filter~ Preda~ 1.21e-2 0.025  0.01206  *        Wilco~
## 2 Lachnospi~ Abunda~ Filter~ Scrap~ 9.73e-4 0.0046 0.00097  ***      Wilco~
## 3 Lachnospi~ Abunda~ Filter~ Shred~ 6.37e-1 0.7    0.63660  ns       Wilco~
## 4 Lachnospi~ Abunda~ Predat~ Scrap~ 2.76e-1 0.34   0.27630  ns       Wilco~
## 5 Lachnospi~ Abunda~ Predat~ Shred~ 2.59e-3 0.0076 0.00259  **       Wilco~
## 6 Lachnospi~ Abunda~ Scraper Shred~ 2.99e-4 0.0042 0.00030  ***      Wilco~
## Filterer Predator  Scraper Shredder 
##      "a"      "b"      "b"      "a" 
## # A tibble: 6 x 9
##   Family     .y.     group1  group2       p  p.adj p.format p.signif method
##   <fct>      <chr>   <chr>   <chr>    <dbl>  <dbl> <chr>    <chr>    <chr> 
## 1 Methyloph~ Abunda~ Filter~ Preda~ 9.15e-1 0.94   0.91511  ns       Wilco~
## 2 Methyloph~ Abunda~ Filter~ Scrap~ 2.08e-3 0.0073 0.00208  **       Wilco~
## 3 Methyloph~ Abunda~ Filter~ Shred~ 1.56e-1 0.22   0.15638  ns       Wilco~
## 4 Methyloph~ Abunda~ Predat~ Scrap~ 8.78e-4 0.0046 0.00088  ***      Wilco~
## 5 Methyloph~ Abunda~ Predat~ Shred~ 2.25e-1 0.3    0.22464  ns       Wilco~
## 6 Methyloph~ Abunda~ Scraper Shred~ 4.57e-4 0.0043 0.00046  ***      Wilco~
## Filterer Predator  Scraper Shredder 
##      "a"      "a"      "b"      "a" 
## # A tibble: 6 x 9
##   Family      .y.     group1  group2       p p.adj p.format p.signif method
##   <fct>       <chr>   <chr>   <chr>    <dbl> <dbl> <chr>    <chr>    <chr> 
## 1 Rhodobacte~ Abunda~ Filter~ Predat~ 0.110  0.17  0.10982  ns       Wilco~
## 2 Rhodobacte~ Abunda~ Filter~ Scraper 0.589  0.67  0.58915  ns       Wilco~
## 3 Rhodobacte~ Abunda~ Filter~ Shredd~ 0.219  0.3   0.21930  ns       Wilco~
## 4 Rhodobacte~ Abunda~ Predat~ Scraper 0.0292 0.051 0.02924  *        Wilco~
## 5 Rhodobacte~ Abunda~ Predat~ Shredd~ 0.0184 0.035 0.01842  *        Wilco~
## 6 Rhodobacte~ Abunda~ Scraper Shredd~ 1      1     1.00000  ns       Wilco~
## Filterer Predator  Scraper Shredder 
##     "ab"      "a"     "ab"      "b" 
## # A tibble: 5 x 9
##   Family    .y.     group1  group2        p  p.adj p.format p.signif method
##   <fct>     <chr>   <chr>   <chr>     <dbl>  <dbl> <chr>    <chr>    <chr> 
## 1 Rikenell~ Abunda~ Filter~ Predat~ 2.58e-2 0.047  0.02577  *        Wilco~
## 2 Rikenell~ Abunda~ Filter~ Scraper 6.67e-3 0.016  0.00667  **       Wilco~
## 3 Rikenell~ Abunda~ Filter~ Shredd~ 3.95e-1 0.46   0.39509  ns       Wilco~
## 4 Rikenell~ Abunda~ Predat~ Shredd~ 2.07e-3 0.0073 0.00207  **       Wilco~
## 5 Rikenell~ Abunda~ Scraper Shredd~ 2.99e-4 0.0042 0.00030  ***      Wilco~
## Filterer Predator  Scraper Shredder 
##      "a"      "b"      "b"      "a" 
## # A tibble: 6 x 9
##   Family     .y.     group1  group2       p  p.adj p.format p.signif method
##   <fct>      <chr>   <chr>   <chr>    <dbl>  <dbl> <chr>    <chr>    <chr> 
## 1 Ruminococ~ Abunda~ Filter~ Preda~ 8.91e-3 0.02   0.00891  **       Wilco~
## 2 Ruminococ~ Abunda~ Filter~ Scrap~ 9.73e-4 0.0046 0.00097  ***      Wilco~
## 3 Ruminococ~ Abunda~ Filter~ Shred~ 1.56e-1 0.22   0.15638  ns       Wilco~
## 4 Ruminococ~ Abunda~ Predat~ Scrap~ 2.76e-1 0.34   0.27630  ns       Wilco~
## 5 Ruminococ~ Abunda~ Predat~ Shred~ 2.59e-3 0.0076 0.00259  **       Wilco~
## 6 Ruminococ~ Abunda~ Scraper Shred~ 2.99e-4 0.0042 0.00030  ***      Wilco~
## Filterer Predator  Scraper Shredder 
##      "a"      "b"      "b"      "a"

GenusAbuAll

Below is a list of the relative abundance of all taxa at the Genus level by Decomp Stage

GenusSorted2
##                Genus      FFG  N        mean          sd          se
## 131 Caldicoprobacter  Scraper  9 0.000000000 0.000000000 0.000000000
## 101      Bacteroides Filterer  4 0.000000000 0.000000000 0.000000000
## 88     Asticcacaulis Shredder  7 0.000000000 0.000000000 0.000000000
## 58       Anaerospora Predator  6 0.133333333 0.183787317 0.075030858
## 40     Agrobacterium Shredder  7 0.004761905 0.012598816 0.004761905
## 129 Caldicoprobacter Filterer  4 0.000000000 0.000000000 0.000000000
## 133      Caulobacter Filterer  4 0.000000000 0.000000000 0.000000000
## 69     Aquabacterium Filterer  4 0.000000000 0.000000000 0.000000000
## 117   Bradyrhizobium Filterer  4 0.000000000 0.000000000 0.000000000
## 100         Bacillus Shredder  7 0.009523810 0.025197632 0.009523810
## 59       Anaerospora  Scraper  9 0.062962963 0.176733041 0.058911014
## 121   Brevibacterium Filterer  4 0.000000000 0.000000000 0.000000000
## 103      Bacteroides  Scraper  9 0.000000000 0.000000000 0.000000000
## 6              02d06 Predator  6 0.000000000 0.000000000 0.000000000
## 70     Aquabacterium Predator  6 0.027777778 0.068041382 0.027777778
## 53      Anaerofustis Filterer  4 0.000000000 0.000000000 0.000000000
## 96      Azospirillum Shredder  7 0.000000000 0.000000000 0.000000000
## 79       Armatimonas  Scraper  9 0.000000000 0.000000000 0.000000000
## 142    Chlorobaculum Predator  6 0.011111111 0.027216553 0.011111111
## 13     Achromobacter Filterer  4 0.000000000 0.000000000 0.000000000
## 119   Bradyrhizobium  Scraper  9 0.122222222 0.259272486 0.086424162
## 77       Armatimonas Filterer  4 0.008333333 0.016666667 0.008333333
## 128     Burkholderia Shredder  7 0.000000000 0.000000000 0.000000000
## 143    Chlorobaculum  Scraper  9 0.000000000 0.000000000 0.000000000
## 111  Bifidobacterium  Scraper  9 0.003703704 0.011111111 0.003703704
## 125     Burkholderia Filterer  4 0.000000000 0.000000000 0.000000000
## 164             CM44 Shredder  7 0.071428571 0.073102089 0.027629992
## 145  Christensenella Filterer  4 0.033333333 0.066666667 0.033333333
## 158      Clostridium Predator 12 0.000000000 0.000000000 0.000000000
## 43         Alistipes  Scraper  9 0.000000000 0.000000000 0.000000000
## 73        Arenimonas Filterer  4 0.008333333 0.016666667 0.008333333
## 27       Actinomyces  Scraper  9 0.000000000 0.000000000 0.000000000
## 82      Arthrobacter Predator  6 0.033333333 0.081649658 0.033333333
## 47      Anaerococcus  Scraper  9 0.003703704 0.011111111 0.003703704
## 141    Chlorobaculum Filterer  4 0.000000000 0.000000000 0.000000000
## 15     Achromobacter  Scraper  9 0.029629630 0.077180244 0.025726748
## 80       Armatimonas Shredder  7 0.047619048 0.066268653 0.025047197
## 17        Acidovorax Filterer  4 0.000000000 0.000000000 0.000000000
## 11               A17  Scraper  9 0.237037037 0.711111111 0.237037037
## 3      [Clostridium]  Scraper  9 0.040740741 0.122222222 0.040740741
## 98          Bacillus Predator  6 0.000000000 0.000000000 0.000000000
## 157      Clostridium Filterer  8 2.845833333 4.455634475 1.575304676
## 159      Clostridium  Scraper 18 0.001851852 0.007856742 0.001851852
## 94      Azospirillum Predator  6 0.138888889 0.340206909 0.138888889
## 140       Cellvibrio Shredder  7 0.042857143 0.113389342 0.042857143
## 31      Actinoplanes  Scraper  9 0.000000000 0.000000000 0.000000000
## 44         Alistipes Shredder  7 0.319047619 0.264475122 0.099962200
## 81      Arthrobacter Filterer  4 0.000000000 0.000000000 0.000000000
## 20        Acidovorax Shredder  7 0.009523810 0.025197632 0.009523810
## 155   Chthoniobacter  Scraper  9 0.096296296 0.233002411 0.077667470
## 54      Anaerofustis Predator  6 0.000000000 0.000000000 0.000000000
## 75        Arenimonas  Scraper  9 0.000000000 0.000000000 0.000000000
## 30      Actinoplanes Predator  6 0.005555556 0.013608276 0.005555556
## 60       Anaerospora Shredder  7 0.142857143 0.280683218 0.106088285
## 55      Anaerofustis  Scraper  9 0.000000000 0.000000000 0.000000000
## 132 Caldicoprobacter Shredder  7 0.028571429 0.040499526 0.015307382
## 149 Chryseobacterium Filterer  4 0.000000000 0.000000000 0.000000000
## 118   Bradyrhizobium Predator  6 0.005555556 0.013608276 0.005555556
## 99          Bacillus  Scraper  9 0.000000000 0.000000000 0.000000000
## 108     Bdellovibrio Shredder  7 0.028571429 0.052453053 0.019825390
## 39     Agrobacterium  Scraper  9 0.000000000 0.000000000 0.000000000
## 35          Afifella  Scraper  9 0.000000000 0.000000000 0.000000000
## 42         Alistipes Predator  6 0.000000000 0.000000000 0.000000000
## 112  Bifidobacterium Shredder  7 0.000000000 0.000000000 0.000000000
## 85     Asticcacaulis Filterer  4 0.000000000 0.000000000 0.000000000
## 24     Acinetobacter Shredder  7 0.047619048 0.063412649 0.023967728
## 29      Actinoplanes Filterer  4 0.008333333 0.016666667 0.008333333
## 151 Chryseobacterium  Scraper  9 0.000000000 0.000000000 0.000000000
## 153   Chthoniobacter Filterer  4 0.025000000 0.031914237 0.015957118
## 137       Cellvibrio Filterer  4 0.000000000 0.000000000 0.000000000
## 90          Azospira Predator  6 0.000000000 0.000000000 0.000000000
## 144    Chlorobaculum Shredder  7 0.000000000 0.000000000 0.000000000
## 56      Anaerofustis Shredder  7 0.014285714 0.037796447 0.014285714
## 126     Burkholderia Predator  6 0.000000000 0.000000000 0.000000000
## 91          Azospira  Scraper  9 0.007407407 0.022222222 0.007407407
## 49       Anaerofilum Filterer  4 1.025000000 0.412198262 0.206099131
## 116  Brachybacterium Shredder  7 0.000000000 0.000000000 0.000000000
## 154   Chthoniobacter Predator  6 0.000000000 0.000000000 0.000000000
## 41         Alistipes Filterer  4 0.000000000 0.000000000 0.000000000
## 61     Anaerotruncus Filterer  4 0.000000000 0.000000000 0.000000000
## 114  Brachybacterium Predator  6 0.000000000 0.000000000 0.000000000
## 115  Brachybacterium  Scraper  9 0.092592593 0.277777778 0.092592593
## 66       Anaerovorax Predator  6 0.000000000 0.000000000 0.000000000
## 124   Brevibacterium Shredder  7 0.000000000 0.000000000 0.000000000
## 113  Brachybacterium Filterer  4 0.000000000 0.000000000 0.000000000
## 163             CM44  Scraper  9 0.000000000 0.000000000 0.000000000
## 19        Acidovorax  Scraper  9 0.000000000 0.000000000 0.000000000
## 92          Azospira Shredder  7 0.000000000 0.000000000 0.000000000
## 161             CM44 Filterer  4 0.000000000 0.000000000 0.000000000
## 21     Acinetobacter Filterer  4 0.383333333 0.589726867 0.294863434
## 38     Agrobacterium Predator  6 0.000000000 0.000000000 0.000000000
## 62     Anaerotruncus Predator  6 0.000000000 0.000000000 0.000000000
## 120   Bradyrhizobium Shredder  7 0.000000000 0.000000000 0.000000000
## 71     Aquabacterium  Scraper  9 0.000000000 0.000000000 0.000000000
## 105     Bdellovibrio Filterer  4 0.158333333 0.142400062 0.071200031
## 150 Chryseobacterium Predator  6 0.027777778 0.068041382 0.027777778
## 1      [Clostridium] Filterer  4 0.000000000 0.000000000 0.000000000
## 86     Asticcacaulis Predator  6 0.011111111 0.027216553 0.011111111
## 34          Afifella Predator  6 0.016666667 0.040824829 0.016666667
## 46      Anaerococcus Predator  6 0.000000000 0.000000000 0.000000000
## 109  Bifidobacterium Filterer  4 0.000000000 0.000000000 0.000000000
## 162             CM44 Predator  6 0.005555556 0.013608276 0.005555556
## 16     Achromobacter Shredder  7 0.000000000 0.000000000 0.000000000
## 68       Anaerovorax Shredder  7 0.342857143 0.359379241 0.135832586
## 89          Azospira Filterer  4 0.000000000 0.000000000 0.000000000
## 102      Bacteroides Predator  6 0.000000000 0.000000000 0.000000000
## 123   Brevibacterium  Scraper  9 0.129629630 0.388888889 0.129629630
## 4      [Clostridium] Shredder  7 0.000000000 0.000000000 0.000000000
## 32      Actinoplanes Shredder  7 0.000000000 0.000000000 0.000000000
## 48      Anaerococcus Shredder  7 0.000000000 0.000000000 0.000000000
## 33          Afifella Filterer  4 0.000000000 0.000000000 0.000000000
## 134      Caulobacter Predator  6 0.000000000 0.000000000 0.000000000
## 104      Bacteroides Shredder  7 0.004761905 0.012598816 0.004761905
## 18        Acidovorax Predator  6 0.000000000 0.000000000 0.000000000
## 52       Anaerofilum Shredder  7 0.519047619 0.304811504 0.115207919
## 130 Caldicoprobacter Predator  6 0.000000000 0.000000000 0.000000000
## 95      Azospirillum  Scraper  9 0.000000000 0.000000000 0.000000000
## 37     Agrobacterium Filterer  4 0.000000000 0.000000000 0.000000000
## 74        Arenimonas Predator  6 0.050000000 0.078173596 0.031914237
## 146  Christensenella Predator  6 0.000000000 0.000000000 0.000000000
## 156   Chthoniobacter Shredder  7 0.004761905 0.012598816 0.004761905
## 9                A17 Filterer  4 0.000000000 0.000000000 0.000000000
## 14     Achromobacter Predator  6 0.000000000 0.000000000 0.000000000
## 67       Anaerovorax  Scraper  9 0.000000000 0.000000000 0.000000000
## 87     Asticcacaulis  Scraper  9 0.000000000 0.000000000 0.000000000
## 136      Caulobacter Shredder  7 0.047619048 0.125988158 0.047619048
## 160      Clostridium Shredder 14 3.611904762 4.123947307 1.102171279
## 51       Anaerofilum  Scraper  9 0.000000000 0.000000000 0.000000000
## 110  Bifidobacterium Predator  6 0.000000000 0.000000000 0.000000000
## 139       Cellvibrio  Scraper  9 0.000000000 0.000000000 0.000000000
## 2      [Clostridium] Predator  6 0.000000000 0.000000000 0.000000000
## 5              02d06 Filterer  4 0.000000000 0.000000000 0.000000000
## 7              02d06  Scraper  9 0.000000000 0.000000000 0.000000000
## 8              02d06 Shredder  7 0.009523810 0.025197632 0.009523810
## 22     Acinetobacter Predator  6 1.933333333 1.262097020 0.515248951
## 23     Acinetobacter  Scraper  9 0.177777778 0.248886409 0.082962136
## 25       Actinomyces Filterer  4 0.000000000 0.000000000 0.000000000
## 36          Afifella Shredder  7 0.000000000 0.000000000 0.000000000
## 83      Arthrobacter  Scraper  9 0.022222222 0.066666667 0.022222222
## 84      Arthrobacter Shredder  7 0.000000000 0.000000000 0.000000000
## 93      Azospirillum Filterer  4 0.000000000 0.000000000 0.000000000
## 106     Bdellovibrio Predator  6 0.244444444 0.223772371 0.091354688
## 107     Bdellovibrio  Scraper  9 0.114814815 0.233399462 0.077799821
## 127     Burkholderia  Scraper  9 0.029629630 0.088888889 0.029629630
## 147  Christensenella  Scraper  9 0.000000000 0.000000000 0.000000000
## 10               A17 Predator  6 0.000000000 0.000000000 0.000000000
## 12               A17 Shredder  7 0.000000000 0.000000000 0.000000000
## 26       Actinomyces Predator  6 0.000000000 0.000000000 0.000000000
## 28       Actinomyces Shredder  7 0.009523810 0.025197632 0.009523810
## 45      Anaerococcus Filterer  4 0.000000000 0.000000000 0.000000000
## 50       Anaerofilum Predator  6 0.000000000 0.000000000 0.000000000
## 57       Anaerospora Filterer  4 0.033333333 0.047140452 0.023570226
## 63     Anaerotruncus  Scraper  9 0.000000000 0.000000000 0.000000000
## 64     Anaerotruncus Shredder  7 0.042857143 0.065868235 0.024895853
## 65       Anaerovorax Filterer  4 0.000000000 0.000000000 0.000000000
## 72     Aquabacterium Shredder  7 0.000000000 0.000000000 0.000000000
## 76        Arenimonas Shredder  7 0.014285714 0.037796447 0.014285714
## 78       Armatimonas Predator  6 0.144444444 0.162845075 0.066481224
## 97          Bacillus Filterer  4 0.000000000 0.000000000 0.000000000
## 122   Brevibacterium Predator  6 0.038888889 0.095257934 0.038888889
## 135      Caulobacter  Scraper  9 0.000000000 0.000000000 0.000000000
## 138       Cellvibrio Predator  6 0.111111111 0.226732017 0.092562958
## 148  Christensenella Shredder  7 0.000000000 0.000000000 0.000000000
## 152 Chryseobacterium Shredder  7 0.033333333 0.074535599 0.028171808

PCoA

Weighted Unifrac

-Ellipses represent 95% CI for the mean of each group

physeqStats<-subset_samples(physeq, FFG!="Predator")
physeqStats<-subset_samples(physeqStats,FFG!="Filterer")
ord=ordinate(physeq,"PCoA", "wunifrac")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))

ord=ordinate(physeq,"PCoA", "wunifrac")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+facet_wrap(~Sampling_station)#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

ord=ordinate(physeq,"PCoA", "wunifrac")
ordplot=plot_ordination(physeq, ord,"samples", color="Sampling_station")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = Sampling_station))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+facet_wrap(~FFG)#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

ord=ordinate(physeq,"PCoA", "wunifrac")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))

Jaccard

-Ellipses represent 95% CI for the mean of each group

ord=ordinate(physeq,"PCoA", "jaccard")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))

ord=ordinate(physeq,"PCoA", "jaccard")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+facet_wrap(~Sampling_station,scales = "free")#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

ord=ordinate(physeq,"PCoA", "jaccard")
ordplot=plot_ordination(physeq, ord,"samples", color="Sampling_station")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = Sampling_station))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+facet_wrap(~FFG)#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

ord=ordinate(physeq,"PCoA", "jaccard")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))

Bray-curtis

-Ellipses represent 95% CI for the mean of each group

ord=ordinate(physeq,"PCoA", "bray")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))

ord=ordinate(physeq,"PCoA", "bray")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+facet_wrap(~Sampling_station)#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

ord=ordinate(physeq,"PCoA", "bray")
ordplot=plot_ordination(physeq, ord,"samples", color="Sampling_station")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = Sampling_station))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+facet_wrap(~FFG)#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

ord=ordinate(physeq,"PCoA", "bray")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))

PERMANOVAs

Weighted Unifrac

Homogenieity of Multivariate Dispersions

physeqStats<-subset_samples(physeq, FFG!="Predator")
physeqStats<-subset_samples(physeqStats,FFG!="Filterer")
GPdist=phyloseq::distance(physeqStats, "wunifrac")
beta=betadisper(GPdist, sample_data(physeqStats)$FFG)
permutest(beta)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df     Sum Sq    Mean Sq      F N.Perm Pr(>F)    
## Groups     1 0.00067157 0.00067157 14.985    999  0.001 ***
## Residuals 14 0.00062744 0.00004482                         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
boxplot(beta)

FFG

GPdist=phyloseq::distance(physeqStats, "wunifrac")
#MONMDS= ordinate(physeq, "NMDS",GPdist)
adonis(GPdist ~ FFG*Sampling_station, as(sample_data(physeqStats), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG * Sampling_station, data = as(sample_data(physeqStats),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                      Df SumsOfSqs   MeanSqs F.Model      R2 Pr(>F)   
## FFG                   1  0.010177 0.0101767  4.9044 0.23867  0.002 **
## Sampling_station      1  0.003762 0.0037616  1.8128 0.08822  0.082 . 
## FFG:Sampling_station  1  0.003800 0.0037998  1.8312 0.08912  0.074 . 
## Residuals            12  0.024900 0.0020750         0.58399          
## Total                15  0.042638                   1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis(GPdist ~ FFG, as(sample_data(physeqStats), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqStats),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs   MeanSqs F.Model      R2 Pr(>F)    
## FFG        1  0.010177 0.0101767   4.389 0.23867  0.001 ***
## Residuals 14  0.032462 0.0023187         0.76133           
## Total     15  0.042638                   1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Jaccard

Homogenieity of Multivariate Dispersions

GPdist=phyloseq::distance(physeq, "jaccard")
beta=betadisper(GPdist, sample_data(physeq)$FFG)
permutest(beta)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)  
## Groups     3 0.048886 0.016295 3.9989    999  0.019 *
## Residuals 22 0.089649 0.004075                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
boxplot(beta)

print("Beta dispersion location")
## [1] "Beta dispersion location"
beta=betadisper(GPdist, sample_data(physeq)$Sampling_station)
permutest(beta)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df    Sum Sq    Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.0017764 0.00177642 2.0475    999  0.161
## Residuals 24 0.0208227 0.00086761
boxplot(beta)

print("Beta dispersion Shredder v Scraper")
## [1] "Beta dispersion Shredder v Scraper"
physeqShredderVScraper<-subset_samples(physeq, FFG=="Shredder"|FFG=="Scraper")


GPdist=phyloseq::distance(physeqShredderVScraper, "jaccard")
beta=betadisper(GPdist, sample_data(physeqShredderVScraper)$FFG)
permutest(beta)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)  
## Groups     1 0.010007 0.010007 3.9245    999  0.071 .
## Residuals 14 0.035700 0.002550                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
boxplot(beta)

print("Beta dispersion Shredder v Predator")
## [1] "Beta dispersion Shredder v Predator"
physeqShredderVPredator<-subset_samples(physeq, FFG=="Shredder"|FFG=="Predator")


GPdist=phyloseq::distance(physeqShredderVPredator, "jaccard")
beta=betadisper(GPdist, sample_data(physeqShredderVPredator)$FFG)
permutest(beta)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df    Sum Sq   Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.0017264 0.0017264 0.9972    999  0.343
## Residuals 11 0.0190437 0.0017312
boxplot(beta)

print("Beta dispersion Shredder v Filterer")
## [1] "Beta dispersion Shredder v Filterer"
physeqShredderVFilterer<-subset_samples(physeq, FFG=="Shredder"|FFG=="Filterer")


GPdist=phyloseq::distance(physeqShredderVFilterer, "jaccard")
beta=betadisper(GPdist, sample_data(physeqShredderVFilterer)$FFG)
permutest(beta)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.016325 0.0163249 2.8041    999  0.155
## Residuals  9 0.052396 0.0058218
boxplot(beta)

print("Beta dispersion Predator v Scraper")
## [1] "Beta dispersion Predator v Scraper"
physeqPredatorVScraper<-subset_samples(physeq, FFG=="Predator"|FFG=="Scraper")


GPdist=phyloseq::distance(physeqPredatorVScraper, "jaccard")
beta=betadisper(GPdist, sample_data(physeqPredatorVScraper)$FFG)
permutest(beta)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq     F N.Perm Pr(>F)
## Groups     1 0.002683 0.0026825 0.936    999  0.347
## Residuals 13 0.037257 0.0028659
boxplot(beta)

print("Beta dispersion Predator v Filterer")
## [1] "Beta dispersion Predator v Filterer"
physeqPredatorVFilterer<-subset_samples(physeq, FFG=="Predator"|FFG=="Filterer")


GPdist=phyloseq::distance(physeqPredatorVFilterer, "jaccard")
beta=betadisper(GPdist, sample_data(physeqPredatorVFilterer)$FFG)
permutest(beta)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq   Mean Sq      F N.Perm Pr(>F)
## Groups     1 0.025561 0.0255605 3.7911    999  0.113
## Residuals  8 0.053937 0.0067422
boxplot(beta)

print("Beta dispersion Scraper v Filterer")
## [1] "Beta dispersion Scraper v Filterer"
physeqScraperVFilterer<-subset_samples(physeq, FFG=="Scraper"|FFG=="Filterer")


GPdist=phyloseq::distance(physeqScraperVFilterer, "jaccard")
beta=betadisper(GPdist, sample_data(physeqScraperVFilterer)$FFG)
permutest(beta)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)  
## Groups     1 0.047159 0.047159 7.3469    999  0.021 *
## Residuals 11 0.070608 0.006419                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
boxplot(beta)

FFG

GPdist=phyloseq::distance(physeqStats, "jaccard")
#MONMDS= ordinate(physeq, "NMDS",GPdist)
adonis(GPdist ~ FFG*Sampling_station, as(sample_data(physeqStats), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG * Sampling_station, data = as(sample_data(physeqStats),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                      Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## FFG                   1    1.2666 1.26661  3.7600 0.19102  0.001 ***
## Sampling_station      1    0.6445 0.64455  1.9133 0.09721  0.012 *  
## FFG:Sampling_station  1    0.6772 0.67721  2.0103 0.10213  0.013 *  
## Residuals            12    4.0424 0.33687         0.60964           
## Total                15    6.6308                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis(GPdist ~ FFG, as(sample_data(physeqStats), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqStats),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## FFG        1    1.2666 1.26661  3.3058 0.19102  0.001 ***
## Residuals 14    5.3642 0.38315         0.80898           
## Total     15    6.6308                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Bray-curtis

Homogenieity of Multivariate Dispersions

GPdist=phyloseq::distance(physeqStats, "bray")
beta=betadisper(GPdist, sample_data(physeqStats)$FFG)
permutest(beta)
## 
## Permutation test for homogeneity of multivariate dispersions
## Permutation: free
## Number of permutations: 999
## 
## Response: Distances
##           Df   Sum Sq  Mean Sq      F N.Perm Pr(>F)  
## Groups     1 0.025993 0.025993 3.8617    999  0.065 .
## Residuals 14 0.094234 0.006731                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
boxplot(beta)

FFG

GPdist=phyloseq::distance(physeqStats, "bray")
#MONMDS= ordinate(physeq, "NMDS",GPdist)
adonis(GPdist ~ FFG*Sampling_station, as(sample_data(physeqStats), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG * Sampling_station, data = as(sample_data(physeqStats),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                      Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## FFG                   1    1.6644 1.66440  6.4304 0.27136  0.001 ***
## Sampling_station      1    0.6540 0.65397  2.5266 0.10662  0.007 ** 
## FFG:Sampling_station  1    0.7092 0.70923  2.7401 0.11563  0.008 ** 
## Residuals            12    3.1060 0.25884         0.50639           
## Total                15    6.1336                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
adonis(GPdist ~ FFG, as(sample_data(physeqStats), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqStats),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## FFG        1    1.6644 1.66440  5.2138 0.27136  0.001 ***
## Residuals 14    4.4692 0.31923         0.72864           
## Total     15    6.1336                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Jaccard Split by Sampling Station

physeqOstana<-subset_samples(physeq, Sampling_station=="Ostana")
physeqOstanaStats<-subset_samples(physeqOstana,FFG!="Predator")
physeqPDRStats<-subset_samples(physeq, Sampling_station!="Ostana")
#MONMDS= ordinate(physeq, "NMDS",GPdist)

Ostana (Predator not included in stats due to N = 2)

GPdist=phyloseq::distance(physeqOstanaStats, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaStats), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaStats),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## FFG        2    1.8065 0.90327  2.6038 0.36653  0.001 ***
## Residuals  9    3.1222 0.34691         0.63347           
## Total     11    4.9287                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print("Shredder vs Scraper Ostana")
## [1] "Shredder vs Scraper Ostana"
physeqOstanaShredderVScraper<-subset_samples(physeqOstanaStats, FFG=="Shredder"|FFG=="Scraper")
GPdist=phyloseq::distance(physeqOstanaShredderVScraper, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaShredderVScraper), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaShredderVScraper),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)  
## FFG        1   0.86139 0.86139  2.3315 0.27984  0.027 *
## Residuals  6   2.21670 0.36945         0.72016         
## Total      7   3.07809                 1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print("Shredder vs Filterer Ostana")
## [1] "Shredder vs Filterer Ostana"
physeqOstanaShredderVFilterer<-subset_samples(physeqOstanaStats, FFG=="Shredder"|FFG=="Filterer")
GPdist=phyloseq::distance(physeqOstanaShredderVFilterer, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaShredderVFilterer), "data.frame"))
## Set of permutations < 'minperm'. Generating entire set.
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaShredderVFilterer),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)  
## FFG        1   0.85695 0.85695  2.8363 0.36194  0.027 *
## Residuals  5   1.51071 0.30214         0.63806         
## Total      6   2.36766                 1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print("Scraper vs Filterer Ostana")
## [1] "Scraper vs Filterer Ostana"
physeqOstanaScraperVFilterer<-subset_samples(physeqOstanaStats, FFG=="Scraper"|FFG=="Filterer")
GPdist=phyloseq::distance(physeqOstanaScraperVFilterer, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaScraperVFilterer), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaScraperVFilterer),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)  
## FFG        1    0.9765 0.97652  2.7159 0.27953  0.011 *
## Residuals  7    2.5169 0.35956         0.72047         
## Total      8    3.4934                 1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PDR (Includes Predator, Shredder, Scraper)

GPdist=phyloseq::distance(physeqPDRStats, "jaccard")

adonis(GPdist ~ FFG, as(sample_data(physeqPDRStats), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRStats),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model     R2 Pr(>F)    
## FFG        2    1.8490 0.92448  2.8964 0.3916  0.001 ***
## Residuals  9    2.8726 0.31918         0.6084           
## Total     11    4.7216                 1.0000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print("Shredder vs Scraper PDR")
## [1] "Shredder vs Scraper PDR"
physeqPDRShredderVScraper<-subset_samples(physeqPDRStats, FFG=="Shredder"|FFG=="Scraper")
GPdist=phyloseq::distance(physeqPDRShredderVScraper, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRShredderVScraper), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRShredderVScraper),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)  
## FFG        1    1.0689 1.06893  3.5129 0.36928  0.029 *
## Residuals  6    1.8257 0.30428         0.63072         
## Total      7    2.8946                 1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print("Shredder vs Predator PDR")
## [1] "Shredder vs Predator PDR"
physeqPDRShredderVPredator<-subset_samples(physeqPDRStats, FFG=="Shredder"|FFG=="Predator")
GPdist=phyloseq::distance(physeqPDRShredderVPredator, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRShredderVPredator), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRShredderVPredator),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)  
## FFG        1   0.90936 0.90936   2.874 0.32387  0.022 *
## Residuals  6   1.89844 0.31641         0.67613         
## Total      7   2.80779                 1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print("Scraper vs Predator PDR")
## [1] "Scraper vs Predator PDR"
physeqPDRScraperVPredator<-subset_samples(physeqPDRStats, FFG=="Scraper"|FFG=="Predator")
GPdist=phyloseq::distance(physeqPDRScraperVPredator, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRScraperVPredator), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRScraperVPredator),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)  
## FFG        1   0.79514 0.79514  2.3605 0.28234   0.02 *
## Residuals  6   2.02107 0.33684         0.71766         
## Total      7   2.81621                 1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Random Forest

By Sampling Station

set.seed(200)
GPr  = transform_sample_counts(physeq, function(x) x / sum(x) ) #transform samples based on relative abundance
GPr = filter_taxa(GPr, function(x) mean(x) > 1e-5, TRUE)
PhylumAll=tax_glom(GPr, "Phylum")

FamilyAll=tax_glom(GPr,"Family")
GenusAll=tax_glom(GPr,"Genus")

ForestData=GenusAll#Change this one so you dont have to rewrite all variables
predictors=t(otu_table(ForestData))
response <- as.factor(sample_data(ForestData)$Sampling_station)
rf.data <- data.frame(response, predictors)
MozzieForest <- randomForest(response~., data = rf.data, ntree = 1000)
print(MozzieForest)#returns overall Random Forest results
## 
## Call:
##  randomForest(formula = response ~ ., data = rf.data, ntree = 1000) 
##                Type of random forest: classification
##                      Number of trees: 1000
## No. of variables tried at each split: 12
## 
##         OOB estimate of  error rate: 3.85%
## Confusion matrix:
##        Ostana PDR class.error
## Ostana     14   0  0.00000000
## PDR         1  11  0.08333333
imp <- importance(MozzieForest)#all the steps that are imp or imp. are building a dataframe that contains info about the taxa used by the Random Forest testto classify treatment 
imp <- data.frame(predictors = rownames(imp), imp)
imp.sort <- arrange(imp, desc(MeanDecreaseGini))
imp.sort$predictors <- factor(imp.sort$predictors, levels = imp.sort$predictors)
imp.20 <- imp.sort[1:23, ]#22
ggplot(imp.20, aes(x = predictors, y = MeanDecreaseGini)) +
  geom_bar(stat = "identity", fill = "indianred") +
  coord_flip() +
  ggtitle("Most important genera for classifying  samples\n by Sampling Station")#\n in a string tells it to start a new line

#imp.20$MeanDecreaseGini
otunames <- imp.20$predictors
r <- rownames(tax_table(ForestData)) %in% otunames
otunames
##  [1] de46410fbc643aeaed03d6aa3878d71a  X1e994910b44683e96c37da4cd04862a4
##  [3] e1a9430f9c3611fc149a9f8e65bf5259  X8766563e48c605e235f4ea2485b02ee3
##  [5] X768d500dec54c80d5ef769902dadd792 e6cad09c4fd44cdb1b919a77f3d92429 
##  [7] X237f7729fcc16a8a864b826f36f307c3 X7b9668ebf6830386ee80c560f7b275a7
##  [9] a046edb519c4e3cbdf77ada497c4d743  X5e419280ad42e77b964625f286da0f08
## [11] e32f997a31db6624fc4f67c269121d4c  X9eab7a1249300ee6e9757c224201f94d
## [13] b4fda0049757d0552cbf3080550ce23b  X2a64f7117c72e249a9c2e5850fd6eb3e
## [15] X147463ced695e1a1dc7c68b5c64c89ed X592e7630b2e58f0bcb2a3f9cccb07ceb
## [17] X887b10e0bdad8789f41fb3c687174051 b672de3b577f55a4350eabdce0f18904 
## [19] a1be228b7d7fe14508a59d112f71fa1a  X129c4eff55474784a99b0c7309c15c50
## [21] ca23f159cd2907c86ea8d9cf10ae5df9  X0fbb434b42e7241efa6451d8f1b429d0
## [23] ad07e5885874179952e16bf2f7bb57b3 
## 165 Levels: de46410fbc643aeaed03d6aa3878d71a ...
PredictorTable<-kable(tax_table(ForestData)[r, ])#returns a list of the most important predictors for Random Forest Classification
PredictorTable
Kingdom Phylum Class Order Family Genus Species Rank5 Rank6 Rank7 Rank4 Rank2 Rank3 Rank1
b672de3b577f55a4350eabdce0f18904 Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Moraxellaceae Perlucidibaca NA NA NA NA NA NA NA NA
b4fda0049757d0552cbf3080550ce23b Bacteria Proteobacteria Alphaproteobacteria Rhodobacterales Rhodobacteraceae Anaerospora NA NA NA NA NA NA NA NA
e1a9430f9c3611fc149a9f8e65bf5259 Bacteria Proteobacteria Alphaproteobacteria Rhodobacterales Rhodobacteraceae Rhodobacter NA NA NA NA NA NA NA NA
ca23f159cd2907c86ea8d9cf10ae5df9 Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Hyphomicrobiaceae Hyphomicrobium NA NA NA NA NA NA NA NA
a046edb519c4e3cbdf77ada497c4d743 Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Rhodoferax NA NA NA NA NA NA NA NA
e6cad09c4fd44cdb1b919a77f3d92429 Bacteria Bacteroidetes [Saprospirae] [Saprospirales] Chitinophagaceae Sediminibacterium NA NA NA NA NA NA NA NA
e32f997a31db6624fc4f67c269121d4c Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Leadbetterella NA NA NA NA NA NA NA NA
ad07e5885874179952e16bf2f7bb57b3 Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Emticicia NA NA NA NA NA NA NA NA
a1be228b7d7fe14508a59d112f71fa1a Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Rudanella NA NA NA NA NA NA NA NA
de46410fbc643aeaed03d6aa3878d71a Bacteria Verrucomicrobia Verrucomicrobiae Verrucomicrobiales Verrucomicrobiaceae Luteolibacter NA NA NA NA NA NA NA NA
GenusRandomForestSubset = subset_taxa(GenusAll, row.names(tax_table(GenusAll))%in% otunames)
GenusRandomForestSubset
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 10 taxa and 26 samples ]
## sample_data() Sample Data:       [ 26 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 10 taxa by 14 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 10 tips and 9 internal nodes ]
df <- psmelt(GenusRandomForestSubset)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Genus", "Sampling_station"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
Trtdata
##                Genus Sampling_station  N       mean         sd          se
## 1        Anaerospora           Ostana 14 0.16666667 0.24564584 0.065651613
## 2        Anaerospora              PDR 12 0.01388889 0.03320683 0.009585986
## 3          Emticicia           Ostana 14 0.50952381 0.54685247 0.146152469
## 4          Emticicia              PDR 12 1.06666667 1.39544714 0.402830892
## 5     Hyphomicrobium           Ostana 14 0.33571429 0.45750425 0.122273153
## 6     Hyphomicrobium              PDR 12 0.18333333 0.40489430 0.116882916
## 7     Leadbetterella           Ostana 14 0.14523810 0.26718128 0.071407201
## 8     Leadbetterella              PDR 12 4.14722222 7.86069956 2.269188505
## 9      Luteolibacter           Ostana 14 0.41904762 0.48739922 0.130262920
## 10     Luteolibacter              PDR 12 0.28333333 0.45979794 0.132732231
## 11     Perlucidibaca           Ostana 14 0.34047619 0.45726399 0.122208941
## 12     Perlucidibaca              PDR 12 1.93055556 3.20673087 0.925703467
## 13       Rhodobacter           Ostana 14 4.97142857 5.19945522 1.389612859
## 14       Rhodobacter              PDR 12 3.25277778 3.30687861 0.954613627
## 15        Rhodoferax           Ostana 14 1.14761905 1.12376917 0.300339943
## 16        Rhodoferax              PDR 12 1.31666667 1.24182173 0.358483055
## 17         Rudanella           Ostana 14 0.16428571 0.36408515 0.097305850
## 18         Rudanella              PDR 12 0.01944444 0.06735753 0.019444444
## 19 Sediminibacterium           Ostana 14 0.05476190 0.08534209 0.022808633
## 20 Sediminibacterium              PDR 12 0.26944444 0.45736583 0.132030142
cdataplot=ggplot(Trtdata, aes(x=Sampling_station,y=mean))+geom_bar(aes(fill = Sampling_station),colour="black", stat="identity")+ facet_grid(~Genus)+xlab("Sampling Station")+ylab("Relative Abundance") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Genus,scales = "free_y")+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
cdataplot

compare_means(Abundance ~ Sampling_station, data = df, group.by = "Genus", p.adjust.method = "fdr",method="kruskal.test")
## # A tibble: 10 x 7
##    Genus             .y.            p p.adj p.format p.signif method       
##    <fct>             <chr>      <dbl> <dbl> <chr>    <chr>    <chr>        
##  1 Leadbetterella    Abundance 0.250   0.5  0.250    ns       Kruskal-Wall~
##  2 Rhodobacter       Abundance 0.136   0.34 0.136    ns       Kruskal-Wall~
##  3 Perlucidibaca     Abundance 0.660   0.73 0.660    ns       Kruskal-Wall~
##  4 Rhodoferax        Abundance 0.796   0.8  0.796    ns       Kruskal-Wall~
##  5 Emticicia         Abundance 0.531   0.66 0.531    ns       Kruskal-Wall~
##  6 Luteolibacter     Abundance 0.0435  0.15 0.044    *        Kruskal-Wall~
##  7 Sediminibacterium Abundance 0.304   0.51 0.304    ns       Kruskal-Wall~
##  8 Hyphomicrobium    Abundance 0.411   0.59 0.411    ns       Kruskal-Wall~
##  9 Rudanella         Abundance 0.0350  0.15 0.035    *        Kruskal-Wall~
## 10 Anaerospora       Abundance 0.0276  0.15 0.028    *        Kruskal-Wall~

ByFFGOstana

GenusAllOstanaStats<-subset_samples(GenusAll, Sampling_station=="Ostana")
GenusAllOstanaStats<-subset_samples(GenusAllOstanaStats, FFG!="Predator")

ForestData=GenusAllOstanaStats#Change this one so you dont have to rewrite all variables
predictors=t(otu_table(ForestData))
response <- as.factor(sample_data(ForestData)$FFG)
rf.data <- data.frame(response, predictors)
MozzieForest <- randomForest(response~., data = rf.data, ntree = 1000)
print(MozzieForest)#returns overall Random Forest results
## 
## Call:
##  randomForest(formula = response ~ ., data = rf.data, ntree = 1000) 
##                Type of random forest: classification
##                      Number of trees: 1000
## No. of variables tried at each split: 12
## 
##         OOB estimate of  error rate: 8.33%
## Confusion matrix:
##          Filterer Scraper Shredder class.error
## Filterer        3       1        0        0.25
## Scraper         0       5        0        0.00
## Shredder        0       0        3        0.00
imp <- importance(MozzieForest)#all the steps that are imp or imp. are building a dataframe that contains info about the taxa used by the Random Forest testto classify treatment 
imp <- data.frame(predictors = rownames(imp), imp)
imp.sort <- arrange(imp, desc(MeanDecreaseGini))
imp.sort$predictors <- factor(imp.sort$predictors, levels = imp.sort$predictors)
imp.20 <- imp.sort[1:19, ]#22
ggplot(imp.20, aes(x = predictors, y = MeanDecreaseGini)) +
  geom_bar(stat = "identity", fill = "indianred") +
  coord_flip() +
  ggtitle("Most important genera for classifying  samples\n by FFG")#\n in a string tells it to start a new line

#imp.20$MeanDecreaseGini
otunames <- imp.20$predictors
r <- rownames(tax_table(ForestData)) %in% otunames
otunames
##  [1] acccb7cec4d146864bc11d37da55dcd0  a7c725b951d62e6f6814fb8ca64a356e 
##  [3] ad07e5885874179952e16bf2f7bb57b3  X7403bfa23d688a94b469051d9f13ae5f
##  [5] X0fbb434b42e7241efa6451d8f1b429d0 b672de3b577f55a4350eabdce0f18904 
##  [7] a42594aa39e8c4b1efb497c275e791ba  a046edb519c4e3cbdf77ada497c4d743 
##  [9] X7b9668ebf6830386ee80c560f7b275a7 dc7d0d97c9604c01c99c972b878e09a0 
## [11] X8766563e48c605e235f4ea2485b02ee3 X07705f7fffe33ad226eab098040fbb75
## [13] X1036caae06a5e5c605ab8120e5b5b7bc e32f997a31db6624fc4f67c269121d4c 
## [15] X43dc593942c5c838cdcb10dac0a39474 d9a11b1b600813533e2f5f0cd7a2eb44 
## [17] X4367180b0dc40f77a063355ebab6a4bb ca23f159cd2907c86ea8d9cf10ae5df9 
## [19] X855eedb949870fe66f915b9293e0de50
## 165 Levels: acccb7cec4d146864bc11d37da55dcd0 ...
PredictorTable<-kable(tax_table(ForestData)[r, ])#returns a list of the most important predictors for Random Forest Classification
PredictorTable
Kingdom Phylum Class Order Family Genus Species Rank5 Rank6 Rank7 Rank4 Rank2 Rank3 Rank1
a7c725b951d62e6f6814fb8ca64a356e Bacteria Firmicutes Clostridia Clostridiales Lachnospiraceae Clostridium NA NA NA NA NA NA NA NA
d9a11b1b600813533e2f5f0cd7a2eb44 Bacteria Firmicutes Clostridia Clostridiales Lachnospiraceae Coprococcus NA NA NA NA NA NA NA NA
dc7d0d97c9604c01c99c972b878e09a0 Bacteria Firmicutes Clostridia Clostridiales Ruminococcaceae Ruminococcus NA NA NA NA NA NA NA NA
b672de3b577f55a4350eabdce0f18904 Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Moraxellaceae Perlucidibaca NA NA NA NA NA NA NA NA
ca23f159cd2907c86ea8d9cf10ae5df9 Bacteria Proteobacteria Alphaproteobacteria Rhizobiales Hyphomicrobiaceae Hyphomicrobium NA NA NA NA NA NA NA NA
a046edb519c4e3cbdf77ada497c4d743 Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Rhodoferax NA NA NA NA NA NA NA NA
acccb7cec4d146864bc11d37da55dcd0 Bacteria Proteobacteria Betaproteobacteria Methylophilales Methylophilaceae Methylotenera NA NA NA NA NA NA NA NA
a42594aa39e8c4b1efb497c275e791ba Bacteria Deferribacteres Deferribacteres Deferribacterales Deferribacteraceae Mucispirillum NA NA NA NA NA NA NA NA
e32f997a31db6624fc4f67c269121d4c Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Leadbetterella NA NA NA NA NA NA NA NA
ad07e5885874179952e16bf2f7bb57b3 Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Emticicia NA NA NA NA NA NA NA NA
GenusRandomForestSubset = subset_taxa(GenusAll, row.names(tax_table(GenusAll))%in% otunames)
GenusRandomForestSubset
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 10 taxa and 26 samples ]
## sample_data() Sample Data:       [ 26 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 10 taxa by 14 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 10 tips and 9 internal nodes ]
df <- psmelt(GenusRandomForestSubset)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Genus", "FFG"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Genus)+xlab("Feeding Group")+ylab("Relative Abundance") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Genus,scales = "free_y")+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
cdataplot

compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr",method="kruskal.test")
## # A tibble: 10 x 7
##    Genus         .y.               p   p.adj p.format p.signif method      
##    <fct>         <chr>         <dbl>   <dbl> <chr>    <chr>    <chr>       
##  1 Leadbetterel~ Abundance 0.000133  0.00027 0.00013  ***      Kruskal-Wal~
##  2 Clostridium   Abundance 0.0000465 0.00023 4.6e-05  ****     Kruskal-Wal~
##  3 Perlucidibaca Abundance 0.0000695 0.00023 7.0e-05  ****     Kruskal-Wal~
##  4 Methylotenera Abundance 0.000339  0.00048 0.00034  ***      Kruskal-Wal~
##  5 Mucispirillum Abundance 0.000268  0.00045 0.00027  ***      Kruskal-Wal~
##  6 Rhodoferax    Abundance 0.000816  0.001   0.00082  ***      Kruskal-Wal~
##  7 Emticicia     Abundance 0.000131  0.00027 0.00013  ***      Kruskal-Wal~
##  8 Ruminococcus  Abundance 0.0000225 0.00023 2.3e-05  ****     Kruskal-Wal~
##  9 Coprococcus   Abundance 0.00663   0.0074  0.00663  **       Kruskal-Wal~
## 10 Hyphomicrobi~ Abundance 0.140     0.14    0.13988  ns       Kruskal-Wal~
Means<-compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr")
Means
## # A tibble: 53 x 9
##    Genus     .y.     group1  group2       p  p.adj p.format p.signif method
##    <fct>     <chr>   <chr>   <chr>    <dbl>  <dbl> <chr>    <chr>    <chr> 
##  1 Leadbett~ Abunda~ Filter~ Preda~ 1.39e-2 0.025  0.01392  *        Wilco~
##  2 Leadbett~ Abunda~ Filter~ Scrap~ 6.53e-3 0.015  0.00653  **       Wilco~
##  3 Leadbett~ Abunda~ Filter~ Shred~ 4.43e-1 0.51   0.44341  ns       Wilco~
##  4 Leadbett~ Abunda~ Predat~ Scrap~ 4.26e-4 0.0047 0.00043  ***      Wilco~
##  5 Leadbett~ Abunda~ Predat~ Shred~ 5.28e-3 0.013  0.00528  **       Wilco~
##  6 Leadbett~ Abunda~ Scraper Shred~ 1.25e-3 0.0052 0.00125  **       Wilco~
##  7 Clostrid~ Abunda~ Filter~ Preda~ 5.74e-3 0.014  0.00574  **       Wilco~
##  8 Clostrid~ Abunda~ Filter~ Scrap~ 9.73e-4 0.0047 0.00097  ***      Wilco~
##  9 Clostrid~ Abunda~ Filter~ Shred~ 7.77e-1 0.79   0.77681  ns       Wilco~
## 10 Clostrid~ Abunda~ Predat~ Shred~ 2.07e-3 0.0061 0.00207  **       Wilco~
## # ... with 43 more rows
Means=compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr")
SigList<-length(unique(Trtdata$Genus))
for (i in levels(Means$Genus)){
  Tax<-i
  TaxAbundance<-subset(Means,Genus==i )
  Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
  difference<-TaxAbundance$p.adj
  names(difference)<-Hyphenated
  Letters<-multcompLetters(difference)
  #print(Letters)
  SigList[i]<-Letters
  
}
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
vec<-unlist(SigList)
vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_wrap(~Genus,scales="free_y")+xlab("FFG")+ylab("Relative Abundance (%) Forest") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+geom_text(aes(x=FFG, y=mean+se+se,label=vec))+ scale_fill_manual(values=cbPalette)
#scale_x_discrete(labels=c("0 hrs", "24 hrs", "48 hrs","72 hrs"))+ theme(axis.text.x = element_text(angle = 45, hjust = 1))#+ scale_fill_manual(values=cbPalette)
cdataplot

ByFFGPDR

GenusAllPDRStats<-subset_samples(GenusAll, Sampling_station!="Ostana")


ForestData=GenusAllPDRStats#Change this one so you dont have to rewrite all variables
predictors=t(otu_table(ForestData))
response <- as.factor(sample_data(ForestData)$FFG)
rf.data <- data.frame(response, predictors)
MozzieForest <- randomForest(response~., data = rf.data, ntree = 1000)
print(MozzieForest)#returns overall Random Forest results
## 
## Call:
##  randomForest(formula = response ~ ., data = rf.data, ntree = 1000) 
##                Type of random forest: classification
##                      Number of trees: 1000
## No. of variables tried at each split: 12
## 
##         OOB estimate of  error rate: 0%
## Confusion matrix:
##          Predator Scraper Shredder class.error
## Predator        4       0        0           0
## Scraper         0       4        0           0
## Shredder        0       0        4           0
imp <- importance(MozzieForest)#all the steps that are imp or imp. are building a dataframe that contains info about the taxa used by the Random Forest testto classify treatment 
imp <- data.frame(predictors = rownames(imp), imp)
imp.sort <- arrange(imp, desc(MeanDecreaseGini))
imp.sort$predictors <- factor(imp.sort$predictors, levels = imp.sort$predictors)
imp.20 <- imp.sort[1:19, ]#22
ggplot(imp.20, aes(x = predictors, y = MeanDecreaseGini)) +
  geom_bar(stat = "identity", fill = "indianred") +
  coord_flip() +
  ggtitle("Most important genera for classifying  samples\n by FFG")#\n in a string tells it to start a new line

#imp.20$MeanDecreaseGini
otunames <- imp.20$predictors
r <- rownames(tax_table(ForestData)) %in% otunames
otunames
##  [1] a046edb519c4e3cbdf77ada497c4d743  X2b05a6e9e6c580f7fc168a5bbb29de2d
##  [3] acccb7cec4d146864bc11d37da55dcd0  ad07e5885874179952e16bf2f7bb57b3 
##  [5] X237f7729fcc16a8a864b826f36f307c3 X5e419280ad42e77b964625f286da0f08
##  [7] X1036caae06a5e5c605ab8120e5b5b7bc X1e994910b44683e96c37da4cd04862a4
##  [9] X07705f7fffe33ad226eab098040fbb75 X6b78b6e57cad2b7f9420cee2c46db2aa
## [11] edb114398ced5bf958f4f404493a6642  e32f997a31db6624fc4f67c269121d4c 
## [13] a7c725b951d62e6f6814fb8ca64a356e  X4367180b0dc40f77a063355ebab6a4bb
## [15] e1a9430f9c3611fc149a9f8e65bf5259  de46410fbc643aeaed03d6aa3878d71a 
## [17] b672de3b577f55a4350eabdce0f18904  dc7d0d97c9604c01c99c972b878e09a0 
## [19] X0fbb434b42e7241efa6451d8f1b429d0
## 165 Levels: a046edb519c4e3cbdf77ada497c4d743 ...
PredictorTable<-kable(tax_table(ForestData)[r, ])#returns a list of the most important predictors for Random Forest Classification
PredictorTable
Kingdom Phylum Class Order Family Genus Species Rank5 Rank6 Rank7 Rank4 Rank2 Rank3 Rank1
a7c725b951d62e6f6814fb8ca64a356e Bacteria Firmicutes Clostridia Clostridiales Lachnospiraceae Clostridium NA NA NA NA NA NA NA NA
dc7d0d97c9604c01c99c972b878e09a0 Bacteria Firmicutes Clostridia Clostridiales Ruminococcaceae Ruminococcus NA NA NA NA NA NA NA NA
b672de3b577f55a4350eabdce0f18904 Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Moraxellaceae Perlucidibaca NA NA NA NA NA NA NA NA
e1a9430f9c3611fc149a9f8e65bf5259 Bacteria Proteobacteria Alphaproteobacteria Rhodobacterales Rhodobacteraceae Rhodobacter NA NA NA NA NA NA NA NA
a046edb519c4e3cbdf77ada497c4d743 Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Rhodoferax NA NA NA NA NA NA NA NA
acccb7cec4d146864bc11d37da55dcd0 Bacteria Proteobacteria Betaproteobacteria Methylophilales Methylophilaceae Methylotenera NA NA NA NA NA NA NA NA
edb114398ced5bf958f4f404493a6642 Bacteria Bacteroidetes [Saprospirae] [Saprospirales] Saprospiraceae Haliscomenobacter NA NA NA NA NA NA NA NA
e32f997a31db6624fc4f67c269121d4c Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Leadbetterella NA NA NA NA NA NA NA NA
ad07e5885874179952e16bf2f7bb57b3 Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Emticicia NA NA NA NA NA NA NA NA
de46410fbc643aeaed03d6aa3878d71a Bacteria Verrucomicrobia Verrucomicrobiae Verrucomicrobiales Verrucomicrobiaceae Luteolibacter NA NA NA NA NA NA NA NA
GenusRandomForestSubset = subset_taxa(GenusAll, row.names(tax_table(GenusAll))%in% otunames)
GenusRandomForestSubset
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 10 taxa and 26 samples ]
## sample_data() Sample Data:       [ 26 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 10 taxa by 14 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 10 tips and 9 internal nodes ]
df <- psmelt(GenusRandomForestSubset)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Genus", "FFG"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Genus)+xlab("Feeding Group")+ylab("Relative Abundance") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Genus,scales = "free_y")+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
#cdataplot

compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr",method="kruskal.test")
## # A tibble: 10 x 7
##    Genus           .y.             p   p.adj p.format p.signif method      
##    <fct>           <chr>       <dbl>   <dbl> <chr>    <chr>    <chr>       
##  1 Leadbetterella  Abundan~  1.33e-4 0.00022 0.00013  ***      Kruskal-Wal~
##  2 Rhodobacter     Abundan~  4.63e-2 0.046   0.04634  *        Kruskal-Wal~
##  3 Clostridium     Abundan~  4.65e-5 0.00022 4.6e-05  ****     Kruskal-Wal~
##  4 Perlucidibaca   Abundan~  6.95e-5 0.00022 7.0e-05  ****     Kruskal-Wal~
##  5 Methylotenera   Abundan~  3.39e-4 0.00048 0.00034  ***      Kruskal-Wal~
##  6 Rhodoferax      Abundan~  8.16e-4 0.001   0.00082  ***      Kruskal-Wal~
##  7 Emticicia       Abundan~  1.31e-4 0.00022 0.00013  ***      Kruskal-Wal~
##  8 Ruminococcus    Abundan~  2.25e-5 0.00022 2.3e-05  ****     Kruskal-Wal~
##  9 Luteolibacter   Abundan~  1.09e-2 0.012   0.01091  *        Kruskal-Wal~
## 10 Haliscomenobac~ Abundan~  1.16e-4 0.00022 0.00012  ***      Kruskal-Wal~
Means<-compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr")
Means
## # A tibble: 55 x 9
##    Genus     .y.     group1  group2       p  p.adj p.format p.signif method
##    <fct>     <chr>   <chr>   <chr>    <dbl>  <dbl> <chr>    <chr>    <chr> 
##  1 Leadbett~ Abunda~ Filter~ Preda~ 1.39e-2 0.025  0.01392  *        Wilco~
##  2 Leadbett~ Abunda~ Filter~ Scrap~ 6.53e-3 0.015  0.00653  **       Wilco~
##  3 Leadbett~ Abunda~ Filter~ Shred~ 4.43e-1 0.51   0.44341  ns       Wilco~
##  4 Leadbett~ Abunda~ Predat~ Scrap~ 4.26e-4 0.0042 0.00043  ***      Wilco~
##  5 Leadbett~ Abunda~ Predat~ Shred~ 5.28e-3 0.014  0.00528  **       Wilco~
##  6 Leadbett~ Abunda~ Scraper Shred~ 1.25e-3 0.0057 0.00125  **       Wilco~
##  7 Rhodobac~ Abunda~ Filter~ Preda~ 1.10e-1 0.15   0.10982  ns       Wilco~
##  8 Rhodobac~ Abunda~ Filter~ Scrap~ 5.89e-1 0.65   0.58915  ns       Wilco~
##  9 Rhodobac~ Abunda~ Filter~ Shred~ 2.18e-1 0.28   0.21825  ns       Wilco~
## 10 Rhodobac~ Abunda~ Predat~ Scrap~ 2.92e-2 0.046  0.02924  *        Wilco~
## # ... with 45 more rows
Means=compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr")
SigList<-length(unique(Trtdata$Genus))
for (i in levels(Means$Genus)){
  Tax<-i
  TaxAbundance<-subset(Means,Genus==i )
  Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
  difference<-TaxAbundance$p.adj
  names(difference)<-Hyphenated
  Letters<-multcompLetters(difference)
  #print(Letters)
  SigList[i]<-Letters
  
}
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
vec<-unlist(SigList)
vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_wrap(~Genus,scales="free_y")+xlab("FFG")+ylab("Relative Abundance (%) Alpine Prairie") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+geom_text(aes(x=FFG, y=mean+se+se,label=vec))+ scale_fill_manual(values=cbPalette)
#scale_x_discrete(labels=c("0 hrs", "24 hrs", "48 hrs","72 hrs"))+ theme(axis.text.x = element_text(angle = 45, hjust = 1))#+ scale_fill_manual(values=cbPalette)
cdataplot

ByFFGCombined

#GenusAllStats<-subset_samples(GenusAll, FFG!="Predator")
#GenusAllStats<-subset_samples(GenusAllStats, FFG!="Filterer")
GenusAllStats<-GenusAll

ForestData=GenusAllStats#Change this one so you dont have to rewrite all variables
predictors=t(otu_table(ForestData))
response <- as.factor(sample_data(ForestData)$FFG)
rf.data <- data.frame(response, predictors)
MozzieForest <- randomForest(response~., data = rf.data, ntree = 1000)
print(MozzieForest)#returns overall Random Forest results
## 
## Call:
##  randomForest(formula = response ~ ., data = rf.data, ntree = 1000) 
##                Type of random forest: classification
##                      Number of trees: 1000
## No. of variables tried at each split: 12
## 
##         OOB estimate of  error rate: 3.85%
## Confusion matrix:
##          Filterer Predator Scraper Shredder class.error
## Filterer        3        1       0        0        0.25
## Predator        0        6       0        0        0.00
## Scraper         0        0       9        0        0.00
## Shredder        0        0       0        7        0.00
imp <- importance(MozzieForest)#all the steps that are imp or imp. are building a dataframe that contains info about the taxa used by the Random Forest testto classify treatment 
imp <- data.frame(predictors = rownames(imp), imp)
imp.sort <- arrange(imp, desc(MeanDecreaseGini))
imp.sort$predictors <- factor(imp.sort$predictors, levels = imp.sort$predictors)
imp.20 <- imp.sort[1:19, ]#22
ggplot(imp.20, aes(x = predictors, y = MeanDecreaseGini)) +
  geom_bar(stat = "identity", fill = "indianred") +
  coord_flip() +
  ggtitle("Most important genera for classifying  samples\n by FFG")#\n in a string tells it to start a new line

#imp.20$MeanDecreaseGini
otunames <- imp.20$predictors
r <- rownames(tax_table(ForestData)) %in% otunames
otunames
##  [1] ad07e5885874179952e16bf2f7bb57b3  X0fbb434b42e7241efa6451d8f1b429d0
##  [3] b672de3b577f55a4350eabdce0f18904  X43dc593942c5c838cdcb10dac0a39474
##  [5] e32f997a31db6624fc4f67c269121d4c  acccb7cec4d146864bc11d37da55dcd0 
##  [7] X4367180b0dc40f77a063355ebab6a4bb a7c725b951d62e6f6814fb8ca64a356e 
##  [9] a046edb519c4e3cbdf77ada497c4d743  dc7d0d97c9604c01c99c972b878e09a0 
## [11] X2b05a6e9e6c580f7fc168a5bbb29de2d a42594aa39e8c4b1efb497c275e791ba 
## [13] edb114398ced5bf958f4f404493a6642  X1036caae06a5e5c605ab8120e5b5b7bc
## [15] X7403bfa23d688a94b469051d9f13ae5f X6b78b6e57cad2b7f9420cee2c46db2aa
## [17] ef53cef3b90f7b7b0189b52eaf440bfd  X237f7729fcc16a8a864b826f36f307c3
## [19] X07705f7fffe33ad226eab098040fbb75
## 165 Levels: ad07e5885874179952e16bf2f7bb57b3 ...
PredictorTable<-kable(tax_table(ForestData)[r, ])#returns a list of the most important predictors for Random Forest Classification
PredictorTable
Kingdom Phylum Class Order Family Genus Species Rank5 Rank6 Rank7 Rank4 Rank2 Rank3 Rank1
a7c725b951d62e6f6814fb8ca64a356e Bacteria Firmicutes Clostridia Clostridiales Lachnospiraceae Clostridium NA NA NA NA NA NA NA NA
dc7d0d97c9604c01c99c972b878e09a0 Bacteria Firmicutes Clostridia Clostridiales Ruminococcaceae Ruminococcus NA NA NA NA NA NA NA NA
b672de3b577f55a4350eabdce0f18904 Bacteria Proteobacteria Gammaproteobacteria Pseudomonadales Moraxellaceae Perlucidibaca NA NA NA NA NA NA NA NA
a046edb519c4e3cbdf77ada497c4d743 Bacteria Proteobacteria Betaproteobacteria Burkholderiales Comamonadaceae Rhodoferax NA NA NA NA NA NA NA NA
acccb7cec4d146864bc11d37da55dcd0 Bacteria Proteobacteria Betaproteobacteria Methylophilales Methylophilaceae Methylotenera NA NA NA NA NA NA NA NA
a42594aa39e8c4b1efb497c275e791ba Bacteria Deferribacteres Deferribacteres Deferribacterales Deferribacteraceae Mucispirillum NA NA NA NA NA NA NA NA
ef53cef3b90f7b7b0189b52eaf440bfd Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Runella NA NA NA NA NA NA NA NA
edb114398ced5bf958f4f404493a6642 Bacteria Bacteroidetes [Saprospirae] [Saprospirales] Saprospiraceae Haliscomenobacter NA NA NA NA NA NA NA NA
e32f997a31db6624fc4f67c269121d4c Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Leadbetterella NA NA NA NA NA NA NA NA
ad07e5885874179952e16bf2f7bb57b3 Bacteria Bacteroidetes Cytophagia Cytophagales Cytophagaceae Emticicia NA NA NA NA NA NA NA NA
GenusRandomForestSubset = subset_taxa(GenusAll, row.names(tax_table(GenusAll))%in% otunames)
GenusRandomForestSubset
## phyloseq-class experiment-level object
## otu_table()   OTU Table:         [ 10 taxa and 26 samples ]
## sample_data() Sample Data:       [ 26 samples by 8 sample variables ]
## tax_table()   Taxonomy Table:    [ 10 taxa by 14 taxonomic ranks ]
## phy_tree()    Phylogenetic Tree: [ 10 tips and 9 internal nodes ]
df <- psmelt(GenusRandomForestSubset)
df$Abundance=df$Abundance*100
Trtdata <- ddply(df, c("Genus", "FFG"), summarise,
                 N    = length(Abundance),
                 mean = mean(Abundance),
                 sd   = sd(Abundance),
                 se   = sd / sqrt(N)
)
Trtdata
##                Genus      FFG N       mean         sd         se
## 1        Clostridium Filterer 4 5.68333333 4.98520032 2.49260016
## 2        Clostridium Predator 6 0.00000000 0.00000000 0.00000000
## 3        Clostridium  Scraper 9 0.00000000 0.00000000 0.00000000
## 4        Clostridium Shredder 7 6.82380952 3.53619980 1.33655789
## 5          Emticicia Filterer 4 0.76666667 0.36004115 0.18002057
## 6          Emticicia Predator 6 2.06666667 1.42439071 0.58150507
## 7          Emticicia  Scraper 9 0.01111111 0.03333333 0.01111111
## 8          Emticicia Shredder 7 0.62380952 0.47442530 0.17931591
## 9  Haliscomenobacter Filterer 4 0.07500000 0.07391186 0.03695593
## 10 Haliscomenobacter Predator 6 0.13888889 0.12003086 0.04900239
## 11 Haliscomenobacter  Scraper 9 0.00000000 0.00000000 0.00000000
## 12 Haliscomenobacter Shredder 7 0.00000000 0.00000000 0.00000000
## 13    Leadbetterella Filterer 4 0.10000000 0.11547005 0.05773503
## 14    Leadbetterella Predator 6 8.38333333 9.64634070 3.93810210
## 15    Leadbetterella  Scraper 9 0.00000000 0.00000000 0.00000000
## 16    Leadbetterella Shredder 7 0.15714286 0.14104673 0.05331065
## 17     Methylotenera Filterer 4 2.50833333 1.45178689 0.72589345
## 18     Methylotenera Predator 6 3.15555556 3.02571693 1.23524377
## 19     Methylotenera  Scraper 9 0.01111111 0.03333333 0.01111111
## 20     Methylotenera Shredder 7 5.09047619 2.19145768 0.82829315
## 21     Mucispirillum Filterer 4 2.76666667 2.05065482 1.02532741
## 22     Mucispirillum Predator 6 0.01666667 0.04082483 0.01666667
## 23     Mucispirillum  Scraper 9 0.00000000 0.00000000 0.00000000
## 24     Mucispirillum Shredder 7 0.81428571 0.57214902 0.21625200
## 25     Perlucidibaca Filterer 4 0.75833333 0.30107831 0.15053915
## 26     Perlucidibaca Predator 6 4.11111111 3.41107130 1.39256403
## 27     Perlucidibaca  Scraper 9 0.00000000 0.00000000 0.00000000
## 28     Perlucidibaca Shredder 7 0.03333333 0.05091751 0.01924501
## 29        Rhodoferax Filterer 4 2.19166667 1.59591492 0.79795746
## 30        Rhodoferax Predator 6 2.16111111 0.97260627 0.39706485
## 31        Rhodoferax  Scraper 9 0.18518519 0.32236587 0.10745529
## 32        Rhodoferax Shredder 7 1.20952381 0.52200266 0.19729846
## 33      Ruminococcus Filterer 4 0.00000000 0.00000000 0.00000000
## 34      Ruminococcus Predator 6 0.00000000 0.00000000 0.00000000
## 35      Ruminococcus  Scraper 9 0.00000000 0.00000000 0.00000000
## 36      Ruminococcus Shredder 7 1.16666667 1.20308246 0.45472243
## 37           Runella Filterer 4 0.27500000 0.44253060 0.22126530
## 38           Runella Predator 6 0.70555556 0.91394546 0.37311667
## 39           Runella  Scraper 9 0.00000000 0.00000000 0.00000000
## 40           Runella Shredder 7 0.01428571 0.03779645 0.01428571
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_grid(~Genus)+xlab("Feeding Group")+ylab("Relative Abundance") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+facet_wrap(~Genus,scales = "free_y")+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+scale_fill_manual(values=cbPalette)
cdataplot

compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr",method="kruskal.test")
## # A tibble: 10 x 7
##    Genus           .y.             p   p.adj p.format p.signif method      
##    <fct>           <chr>       <dbl>   <dbl> <chr>    <chr>    <chr>       
##  1 Leadbetterella  Abundan~  1.33e-4 0.00022 0.00013  ***      Kruskal-Wal~
##  2 Clostridium     Abundan~  4.65e-5 0.00022 4.6e-05  ****     Kruskal-Wal~
##  3 Perlucidibaca   Abundan~  6.95e-5 0.00022 7.0e-05  ****     Kruskal-Wal~
##  4 Methylotenera   Abundan~  3.39e-4 0.00038 0.00034  ***      Kruskal-Wal~
##  5 Mucispirillum   Abundan~  2.68e-4 0.00036 0.00027  ***      Kruskal-Wal~
##  6 Rhodoferax      Abundan~  8.16e-4 0.00082 0.00082  ***      Kruskal-Wal~
##  7 Emticicia       Abundan~  1.31e-4 0.00022 0.00013  ***      Kruskal-Wal~
##  8 Ruminococcus    Abundan~  2.25e-5 0.00022 2.3e-05  ****     Kruskal-Wal~
##  9 Runella         Abundan~  2.89e-4 0.00036 0.00029  ***      Kruskal-Wal~
## 10 Haliscomenobac~ Abundan~  1.16e-4 0.00022 0.00012  ***      Kruskal-Wal~
Means<-compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr")
Means
## # A tibble: 55 x 9
##    Genus     .y.     group1  group2       p  p.adj p.format p.signif method
##    <fct>     <chr>   <chr>   <chr>    <dbl>  <dbl> <chr>    <chr>    <chr> 
##  1 Leadbett~ Abunda~ Filter~ Preda~ 1.39e-2 0.022  0.01392  *        Wilco~
##  2 Leadbett~ Abunda~ Filter~ Scrap~ 6.53e-3 0.013  0.00653  **       Wilco~
##  3 Leadbett~ Abunda~ Filter~ Shred~ 4.43e-1 0.48   0.44341  ns       Wilco~
##  4 Leadbett~ Abunda~ Predat~ Scrap~ 4.26e-4 0.0036 0.00043  ***      Wilco~
##  5 Leadbett~ Abunda~ Predat~ Shred~ 5.28e-3 0.012  0.00528  **       Wilco~
##  6 Leadbett~ Abunda~ Scraper Shred~ 1.25e-3 0.0046 0.00125  **       Wilco~
##  7 Clostrid~ Abunda~ Filter~ Preda~ 5.74e-3 0.012  0.00574  **       Wilco~
##  8 Clostrid~ Abunda~ Filter~ Scrap~ 9.73e-4 0.0041 0.00097  ***      Wilco~
##  9 Clostrid~ Abunda~ Filter~ Shred~ 7.77e-1 0.79   0.77681  ns       Wilco~
## 10 Clostrid~ Abunda~ Predat~ Shred~ 2.07e-3 0.0054 0.00207  **       Wilco~
## # ... with 45 more rows
Means=compare_means(Abundance ~ FFG, data = df, group.by = "Genus", p.adjust.method = "fdr")
SigList<-length(unique(Trtdata$Genus))
for (i in levels(Means$Genus)){
  Tax<-i
  TaxAbundance<-subset(Means,Genus==i )
  Hyphenated<-as.character(paste0(TaxAbundance$group1,"-",TaxAbundance$group2))
  difference<-TaxAbundance$p.adj
  names(difference)<-Hyphenated
  Letters<-multcompLetters(difference)
  #print(Letters)
  SigList[i]<-Letters
  
}
## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length

## Warning in SigList[i] <- Letters: number of items to replace is not a
## multiple of replacement length
vec<-unlist(SigList)
vec<-vec[-1]
cdataplot=ggplot(Trtdata, aes(x=FFG,y=mean))+geom_bar(aes(fill = FFG),colour="black", stat="identity")+ facet_wrap(~Genus,scales="free_y")+xlab("FFG")+ylab("Relative Abundance (%)") + theme(axis.text.x = element_text(angle = 45, hjust = 1))+theme(axis.title.x=element_blank())+geom_errorbar(aes(ymin=mean-se,ymax=mean+se))+geom_text(aes(x=FFG, y=mean+se+se,label=vec))+ scale_fill_manual(values=cbPalette)
#scale_x_discrete(labels=c("0 hrs", "24 hrs", "48 hrs","72 hrs"))+ theme(axis.text.x = element_text(angle = 45, hjust = 1))#+ scale_fill_manual(values=cbPalette)
cdataplot

PiCRUST

PiCRUST<-import_biom("C:\\Users\\Joe Receveur\\Downloads\\ItalyInvertPiCRUST2.20.19.biom")
PiCRUST<-merge_phyloseq(PiCRUST,sampdat)
PiCRUST=rarefy_even_depth(PiCRUST, 750000, replace = TRUE, trimOTUs = TRUE, verbose = TRUE,rngseed = TRUE)
## `set.seed(TRUE)` was used to initialize repeatable random subsampling.
## Please record this for your records so others can reproduce.
## Try `set.seed(TRUE); .Random.seed` for the full vector
## ...
## 487OTUs were removed because they are no longer 
## present in any sample after random subsampling
## ...

PiCRUST by sampling x FFG

GPdist=phyloseq::distance(PiCRUST, "jaccard")
adonis(GPdist ~ FFG*Sampling_station, as(sample_data(PiCRUST), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG * Sampling_station, data = as(sample_data(PiCRUST),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##                      Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)    
## FFG                   3   0.90035 0.300117  5.0869 0.40838  0.001 ***
## Sampling_station      1   0.08310 0.083101  1.4085 0.03769  0.175    
## FFG:Sampling_station  2   0.10029 0.050145  0.8499 0.04549  0.641    
## Residuals            19   1.12096 0.058998         0.50844           
## Total                25   2.20470                  1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PiCRUST beta div by FFG

print("Shredder vs Scraper")
## [1] "Shredder vs Scraper"
PiCRUSTShredderVScraper<-subset_samples(PiCRUST, FFG=="Shredder"|FFG=="Scraper")
GPdist=phyloseq::distance(PiCRUSTShredderVScraper, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(PiCRUSTShredderVScraper), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(PiCRUSTShredderVScraper),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
## FFG        1   0.50616 0.50616  6.4478 0.31533  0.001 ***
## Residuals 14   1.09902 0.07850         0.68467           
## Total     15   1.60518                 1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print("Shredder vs Filterer")
## [1] "Shredder vs Filterer"
PiCRUSTShredderVFilterer<-subset_samples(PiCRUST, FFG=="Shredder"|FFG=="Filterer")
GPdist=phyloseq::distance(PiCRUSTShredderVFilterer, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(PiCRUSTShredderVFilterer), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(PiCRUSTShredderVFilterer),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs  MeanSqs F.Model     R2 Pr(>F)  
## FFG        1  0.057191 0.057191  2.5935 0.2237  0.027 *
## Residuals  9  0.198462 0.022051         0.7763         
## Total     10  0.255653                  1.0000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print("Shredder vs Predator")
## [1] "Shredder vs Predator"
physeqShredderVPredator<-subset_samples(PiCRUST, FFG=="Shredder"|FFG=="Predator")
GPdist=phyloseq::distance(physeqShredderVPredator, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqShredderVPredator), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqShredderVPredator),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
## FFG        1   0.33225 0.33225  17.331 0.61173  0.002 **
## Residuals 11   0.21088 0.01917         0.38827          
## Total     12   0.54313                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print("Scraper vs Filterer")
## [1] "Scraper vs Filterer"
PiCRUSTScraperVFilterer<-subset_samples(PiCRUST, FFG=="Scraper"|FFG=="Filterer")
GPdist=phyloseq::distance(PiCRUSTScraperVFilterer, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(PiCRUSTScraperVFilterer), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(PiCRUSTScraperVFilterer),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)  
## FFG        1   0.29397 0.293967  2.9572 0.21188  0.014 *
## Residuals 11   1.09347 0.099406         0.78812         
## Total     12   1.38744                  1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print("Scraper vs Predator")
## [1] "Scraper vs Predator"
physeqScraperVPredator<-subset_samples(PiCRUST, FFG=="Scraper"|FFG=="Predator")
GPdist=phyloseq::distance(physeqScraperVPredator, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqScraperVPredator), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqScraperVPredator),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)   
## FFG        1    0.3455 0.34550  4.0614 0.23805  0.002 **
## Residuals 13    1.1059 0.08507         0.76195          
## Total     14    1.4514                 1.00000          
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print("Predator vs Filterer")
## [1] "Predator vs Filterer"
PiCRUSTPredatorVFilterer<-subset_samples(PiCRUST, FFG=="Predator"|FFG=="Filterer")
GPdist=phyloseq::distance(PiCRUSTPredatorVFilterer, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(PiCRUSTPredatorVFilterer), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(PiCRUSTPredatorVFilterer),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)  
## FFG        1   0.13581 0.135811  5.2914 0.39811  0.013 *
## Residuals  8   0.20533 0.025666         0.60189         
## Total      9   0.34114                  1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PiCRUST Jaccard Split by Sampling Station

physeqOstana<-subset_samples(PiCRUST, Sampling_station=="Ostana")
physeqOstanaStats<-subset_samples(physeqOstana,FFG!="Predator")
physeqPDRStats<-subset_samples(PiCRUST, Sampling_station!="Ostana")
#MONMDS= ordinate(physeq, "NMDS",GPdist)

Ostana (Predator not included in stats due to N = 2)

GPdist=phyloseq::distance(physeqOstanaStats, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaStats), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaStats),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs  MeanSqs F.Model     R2 Pr(>F)  
## FFG        2   0.32232 0.161160  2.2914 0.3374  0.018 *
## Residuals  9   0.63299 0.070333         0.6626         
## Total     11   0.95531                  1.0000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print("Shredder vs Scraper Ostana")
## [1] "Shredder vs Scraper Ostana"
physeqOstanaShredderVScraper<-subset_samples(physeqOstanaStats, FFG=="Shredder"|FFG=="Scraper")
GPdist=phyloseq::distance(physeqOstanaShredderVScraper, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaShredderVScraper), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaShredderVScraper),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)
## FFG        1   0.20183 0.201826   2.257 0.27334   0.11
## Residuals  6   0.53654 0.089423         0.72666       
## Total      7   0.73836                  1.00000
print("Shredder vs Filterer Ostana")
## [1] "Shredder vs Filterer Ostana"
physeqOstanaShredderVFilterer<-subset_samples(physeqOstanaStats, FFG=="Shredder"|FFG=="Filterer")
GPdist=phyloseq::distance(physeqOstanaShredderVFilterer, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaShredderVFilterer), "data.frame"))
## Set of permutations < 'minperm'. Generating entire set.
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaShredderVFilterer),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)
## FFG        1   0.03489 0.034890  1.2951 0.20573  0.294
## Residuals  5   0.13470 0.026939         0.79427       
## Total      6   0.16959                  1.00000
print("Scraper vs Filterer Ostana")
## [1] "Scraper vs Filterer Ostana"
physeqOstanaScraperVFilterer<-subset_samples(physeqOstanaStats, FFG=="Scraper"|FFG=="Filterer")
GPdist=phyloseq::distance(physeqOstanaScraperVFilterer, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqOstanaScraperVFilterer), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqOstanaScraperVFilterer),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs  MeanSqs F.Model     R2 Pr(>F)  
## FFG        1   0.22322 0.223224  2.6273 0.2729  0.024 *
## Residuals  7   0.59475 0.084965         0.7271         
## Total      8   0.81798                  1.0000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PDR (Includes Predator, Shredder, Scraper)

GPdist=phyloseq::distance(physeqPDRStats, "jaccard")

adonis(GPdist ~ FFG, as(sample_data(physeqPDRStats), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRStats),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)    
## FFG        2   0.61094 0.305470  6.0445 0.57324  0.001 ***
## Residuals  9   0.45483 0.050537         0.42676           
## Total     11   1.06577                  1.00000           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print("Shredder vs Scraper PDR")
## [1] "Shredder vs Scraper PDR"
physeqPDRShredderVScraper<-subset_samples(physeqPDRStats, FFG=="Shredder"|FFG=="Scraper")
GPdist=phyloseq::distance(physeqPDRShredderVScraper, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRShredderVScraper), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRShredderVScraper),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)  
## FFG        1   0.35923 0.35923  5.2268 0.46556  0.041 *
## Residuals  6   0.41237 0.06873         0.53444         
## Total      7   0.77161                 1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print("Shredder vs Predator PDR")
## [1] "Shredder vs Predator PDR"
physeqPDRShredderVPredator<-subset_samples(physeqPDRStats, FFG=="Shredder"|FFG=="Predator")
GPdist=phyloseq::distance(physeqPDRShredderVPredator, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRShredderVPredator), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRShredderVPredator),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)  
## FFG        1   0.27869 0.278694  24.559 0.80366   0.04 *
## Residuals  6   0.06809 0.011348         0.19634         
## Total      7   0.34678                  1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
print("Scraper vs Predator PDR")
## [1] "Scraper vs Predator PDR"
physeqPDRScraperVPredator<-subset_samples(physeqPDRStats, FFG=="Scraper"|FFG=="Predator")
GPdist=phyloseq::distance(physeqPDRScraperVPredator, "jaccard")
adonis(GPdist ~ FFG, as(sample_data(physeqPDRScraperVPredator), "data.frame"))
## 
## Call:
## adonis(formula = GPdist ~ FFG, data = as(sample_data(physeqPDRScraperVPredator),      "data.frame")) 
## 
## Permutation: free
## Number of permutations: 999
## 
## Terms added sequentially (first to last)
## 
##           Df SumsOfSqs  MeanSqs F.Model      R2 Pr(>F)  
## FFG        1   0.27848 0.278483  3.8931 0.39352  0.028 *
## Residuals  6   0.42920 0.071533         0.60648         
## Total      7   0.70768                  1.00000         
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

PCoA PiCRUST

ord=ordinate(PiCRUST,"PCoA", "jaccard")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())#+ theme(legend.justification=c(1,0), legend.position=c(1,0))

ord=ordinate(PiCRUST,"PCoA", "jaccard")
ordplot=plot_ordination(physeq, ord,"samples", color="FFG")+geom_point(size=4)+scale_colour_manual(values=cbPalette)+scale_fill_manual(values=cbPalette)
ordplot+ stat_ellipse(type= "norm",geom = "polygon", alpha = 1/4, aes(fill = FFG))+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+facet_wrap(~Sampling_station)#+ theme(legend.justification=c(1,0), legend.position=c(1,0))
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse